# Analysis of thermal effects in a cavitating inducer using Rayleigh equation

Jean-Pierre Franc, Christian Pellone

<figure><img src="image_1.png" /></figure>

## To cite this version:

Jean-Pierre Franc, Christian Pellone. Analysis of thermal effects in a cavitating inducer using Rayleigh equation. Journal of Fluids Engineering, 2007, 129, pp.974-983. 10.1115/1.2746919. hal-00265809

HAL Id: hal-00265809

https://hal.science/hal-00265809v1

Submitted on 13 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Analysis of Thermal Effects in a Cavitating Inducer Using Rayleigh Equation

A simple model based on the resolution of Rayleigh equation is used to analyze thermal effects in cavitation. Two different assumptions are considered for the modeling of heat transfer toward the liquid/vapor interface. One is based upon a convective type approach using a convection heat transfer coefficient or the equivalent Nusselt number. The other one is based upon the resolution of the heat diffusion equation in the liquid surrounding the bubble. This conductive-type approach requires one to specify the eddy thermal diffusivity or the equivalent Peclet number. Both models are applied to a cavitating inducer. The basic pressure distribution on the blades is determined from a potential flow computation in a two-dimensional cascade of flat plates. The sheet cavity, which develops from the leading edge, is approximated by the envelope of a hemispherical bubble traveling on the suction side of the blade. Cavity shape and temperature distribution predicted by both models are compared. The evolutions of cavity length with the cavitation number for cold water (without thermal effects) and for Refrigerant 114 at two different temperatures is compared to experimental data. Such a simple model is easy to apply and appears to be quite pertinent for the analysis of thermal effects in a cavitating inducer.

## 1 Introduction

In a number of fluids, the development of cavitation goes with significant thermal effects. The so-called thermodynamic effect has been extensively studied since the early works of Stahl et al. [1] and Stepanoff [2]. The phenomenon is characterized, in particular, by a temperature drop inside two-phase regions. It is due to the latent heat of vaporization, which is taken to the liquid surrounding the cavitating zone whose temperature is then lowered.

The effect is negligible in water at room temperature but may be significant at a higher temperature or in other liquids as, e.g., cryogenic fluids used for the propulsion of space rockets. Several parameters have been introduced to estimate the importance of thermal effects in cavitating flows. Stepanoff [2] defined a characteristic temperature drop by

$$ \Delta T^* = \frac{\rho_v L}{\rho_\ell c_{p\ell}} \qquad (1) $$

This is the temperature drop to be applied to a unit volume of liquid to supply the latent heat required for the vaporization of a unit volume of vapor. Experiments show that this parameter is a relevant order of magnitude of the temperature drop in cavitating flows. For water at room temperature, it is quite small and equal to 0.01 K, whereas it exceeds 1 K for liquid hydrogen at 22 K. This big difference demonstrates the major influence of the liquid properties on the phenomenon. The well-known B factor of Stepanoff [2] is the nondimensional parameter:

$$ B = \frac{\Delta T}{\Delta T^*} \qquad (2) $$

where ΔT is the actual temperature drop in the cavitating flow.

Using the slope $dp_v/dT$ of the vapor pressure curve, a characteristic pressure difference can be built on the basis of the previous characteristic temperature difference ΔT*:

Jean-Pierre Franc

LEGI, BP 53,
38041 Grenoble Cedex 9, France

Christian Pellone

LEGI, BP 53,
38041 Grenoble Cedex 9, France

$$ \Delta p_v^* = \frac{dp_v}{dT} \Delta T^* \qquad (3) $$

This is the drop in vapor pressure associated with the drop in temperature ΔT*. Since ΔT* is a relevant order of magnitude of the temperature drop, Δp_v*, can be considered as a relevant order of magnitude of the pressure drop inside the cavities. This is true as long as the thermodynamic equilibrium condition is achieved, which is generally the case. Nonequilibrium is expected only at the very final stage of collapse when the kinetics of phase change can become a limiting factor because of the very high velocity of the interface.

Brennen [3,4] suggests the use of the thermodynamic parameter

$$ \Sigma = \frac{\Delta p_v^*}{\rho_\ell \sqrt{\alpha_\ell}} \qquad (4) $$

to measure the importance of thermal effects in cavitating flows. The Σ parameter was derived from the analysis of the dynamics of a single bubble with the assumption that heat transfer to the interface is provided by conduction through the surrounding liquid. The thickness of the thermal boundary layer then plays an important role and accounts for the presence of the square root of the liquid thermal diffusivity α_ℓ. The Σ parameter has the units m/s³/².

Considering the case of a sheet cavity and still assuming heat transfer to the cavity interface by conduction through the neighboring liquid, Kato [5] introduced a similar parameter defined by

$$ \alpha = \sqrt{\frac{\rho_\ell \Sigma}{\rho_v}} \qquad (5) $$

The importance of thermal effects increases with α and Σ. From a systematic analysis of Hord's experiments [6], Kato [5] showed a good correlation of the nondimensional pressure drop $2Δp_v/ρ_\ell V^2$ inside the cavity with the nondimensional parameter $α\sqrt{c/V^3}$, where c and V are a reference length scale and a reference velocity scale of the flow.

It is of major concern to observe that, in such models based on a well-defined interface between the liquid and either a cavitation bubble or a sheet cavity, heat transfer by conduction characterized
<figure><img src="image_3.png" /><figcaption>Fig. 1 Scheme of the 2D cascade (vertical scale has been expanded)</figcaption></figure>

by the classical liquid thermal diffusivity $\alpha_{\ell}$ is far too small to account for the measured values of the temperature drop. For a correct prediction, Kato [5] and Watanabe et al. [7] have shown that the liquid thermal diffusivity should be multiplied by an am- plifying coefficient $\varepsilon$ of the order of 8000 [7] up to $10^5$ [5]. This large value is supposed to account for the enhancement of heat transfer due to both turbulence and an increase in the interface area probably due to the small-scale instabilities of the interface. The importance of the so-called interface roughness was con- firmed by Fruman et al. [8,9] on the basis of Billet’s entrainment theory [10–12].

Besides interface models, two-phase models have been devel- oped. The two-phase mixture is generally considered as a pseud- ofluid, which obeys the Navier–Stokes equations. A suitable con- stitutive equation for the mixture as the barotropic law can be specified (Rapposelli et d’Agostino [13], Rolland et al. [14]). An- other option consists in solving an additional transport equation for the vapor phase with a source term modeling the cavitation phenomenon (see, e.g., Tani and Nagashima [15] and Hosangadi and Ahuja [16–18]). The vaporization or condensation process can be modeled either by ad hoc equations as in Ref. [16] or by seeding the flow field with microbubbles and computing their evo- lution using the Rayleigh–Plesset equation, as proposed originally by Kubota et al. [19] (see, e.g., Singhal et al. [20]). To address thermal effects, it is necessary to solve the energy equation for the mixture together with the mass and momentum conservation equations. These equations are essentially coupled through the dependence of the vapor pressure with the temperature.

For interface models, the volume of liquid surrounding vapor is almost infinite and heat transfer is supposed to be limited by the growth rate of the thermal boundary layer. Conversely, two-phase models assume an intimate mixing of liquid and vapor and the volume of liquid surrounding a given bubble is naturally limited by the surrounding bubbles. Hence, thermal effects are essentially regulated by the limited volume of liquid able to supply the heat for vaporization and not by heat diffusion as for interface models. Thermal effects are then basically controlled by the evolution of void fraction without any influence of the liquid thermal diffusiv- ity.

The difference between both types of models can easily be demonstrated by comparing basic expressions of the B factor. In the case of a bubble whose radius grows from 0 to R during time t, the B factor is given by (see, e.g., Ref. [21])

$$B \approx \frac{R}{\sqrt{\alpha_{\ell} t}} \qquad (6)$$

On the other hand, for a two-phase mixture of void fraction $\alpha_v$, the heat conservation equation leads to (see Ref. [22]):

$$B \approx \frac{\alpha_v}{1 - \alpha_v} \qquad (7)$$

In the case of Eq. (6), thermal effects depend upon the liquid thermal diffusivity, which is obviously not the case for Eq. (7). In addition, Eq. (7) becomes singular when the void fraction ap- proaches unity, which shows that the two-phase model is valid for rather dilute mixtures but fails when the cavitating region be- comes pure vapor. If so, it is necessary to switch to interface-type models. Hence, it is of primary importance, in practice, to recog- nize the type of cavities in order to be able to identify the appro- priate physical mechanism, which actually limits the supply of heat necessary for vaporization.

The present paper is devoted to an analysis of thermal effects on the basis of the Rayleigh equation. It belongs to the class of interface models for which heat transfer through the thermal boundary layer is supposed to be the limiting physical phenom- enon. The present analysis is very similar to that conducted by Watanabe et al. [7] except that it considers bubble cavitation in- stead of sheet cavitation. It is also applied to the case of a cavi- tating inducer of a space rocket turbopump investigated experi- mentally by Franc et al. [22].

The analysis comprises two successive steps. First, the pressure distribution in a two-dimensional (2D) noncavitating cascade is computed using a boundary element method. Second, the Ray- leigh equation is solved using the previous pressure distribution in order to compute the evolution of a cavitation bubble on the blade suction side.

The original aspect of the present contribution essentially lies in the thermal term added to the Rayleigh equation in order to ac- count for thermal aspects. Two different forms based either on a convective- or a conductive-type approach are considered. Both approaches are evaluated by comparison with the experimental results [22]. The main objective is to evaluate whether this kind of simple model is able to account for the trends observed experi- mentally as the increase of the B factor with the cavity length.

In the present approach, the sheet cavity observed experimen- tally is compared to the computed envelope of a hemispherical bubble, which grows and collapses on the suction side of the blades. This kind of simple model obviously suffers from limita- tions. One is the lack of any feedback of cavitation on the pressure distribution since the bubble evolution is computed using the non- cavitating pressure distribution. The 2D nature of the computation of the pressure distribution is also a limitation since inducers may be subjected to backflow, which may significantly change angles of attack. Nevertheless, the present work supplies a theoretical support for the physical interpretation and understanding of ex- perimental results.
<figure><img src="image_4.png" /><figcaption>Fig. 2 Close view of the blade leading edge</figcaption></figure>

## 2 Pressure Distribution Computation

The inducer is represented by a 2D cascade of flat plates, as shown in Fig. 1. This 2D configuration is obtained from a cylindrical section of the inducer at radius $r$. Enlargements of the leading edge and trailing edge are shown in Figs. 2 and 3. Calculation of the bubble evolution was started at the back end of the bevel, from point B where the pressure is minimum (see also Fig. 4). This location corresponds to the observed point of detachment of cavities. As for the trailing edge, its shape was chosen in order to allow a smooth detachment of the flow in a direction parallel to the blades.

The relative thickness of the blade is about 1% and the solidity $c/H$ is 2. The incidence angle $\alpha_a$ and the stagger angle $\gamma$ of the cascade are such that the blade bevel is parallel to the j direction of the grid. The inlet velocity $V_u$ is supposed perpendicular to the grid (no prerotation). At upstream infinity, the reference relative velocity is

$$W_u = \sqrt{V_u^2 + \omega^2 r^2} \quad (8)$$

where $\omega$ is the angular speed of rotation. Any radial flow normal to the cylinder surface is ignored and the flow is then supposed to be 2D in the (i,j) plane. In addition, it is assumed incompressible, irrotational, and inviscid. Under these assumptions, the pressure coefficient is given by

$$C_p = 1 - \frac{W^2 - \omega^2 r^2}{W_u^2 - \omega^2 r^2} \quad (9)$$

The pressure distribution on the blade is computed using the boundary element method. A double distribution of sources and

<figure><img src="image_5.png" /><figcaption>Fig. 3 Close view of the blade trailing edge</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4 Pressure distribution on the blades. Comparison between the present results and that of Watanabe et al. (Ref. [7]).</figcaption></figure>

vortices ($\sigma^*, \gamma^*$) is distributed on the blade surface with the space periodicity $H$ in the j direction to take into account the infinite series of blades. This distribution is chosen so that the relative velocity at upstream infinity is in the I direction. The boundary conditions consist in a zero normal velocity on the blade (slip condition) and the continuity of pressure at the trailing edge. Conservation of flowrate through the inducer requires that the component of the relative velocity $W \cdot i$ is the same at upstream and downstream infinity. The flow at downstream infinity is deviated with an angle $\phi$ given by

$$\tan \phi = \frac{(\sin \gamma) \oint_{\text{blade}} \gamma^* dx}{HW_u + (\cos \gamma) \oint_{\text{blade}} \gamma^* dx} \quad (10)$$

The double distribution used here introduces more unknowns than equations to be solved. Additional conditions are imposed in order to minimize the discretization errors. This allows us to compute, with an improved accuracy, the pressure coefficient on blades of small relative thickness.

The computational method has been validated by comparison with the results obtained by Watanabe et al. [7]. Both computations agree fairly well, as shown in Fig. 4. The differences, which are observed particularly at the leading edge near the bevel, are due to the nonlinear nature of the present computation contrary to Ref. [7], which assumes blades of zero thickness.

## 3 Basic Equations

### 3.1 Nondimensional Form of the Rayleigh Equation. Once the pressure distribution has been computed, Rayleigh equation

$$R \frac{d^2 R}{dt^2} + \frac{3}{2} \left( \frac{dR}{dt} \right)^2 = \frac{p_v(T_c) - p(t)}{\rho_\ell} \quad (11)$$

is solved in order to determine the evolution of the radius $R$ of a bubble traveling on the suction side of the blades. $p(t)$ is the time dependent pressure to which the bubble is subjected as it moves along the blades. The effects of viscosity and gas content, which generally become negligible as soon as the original microbubble becomes a macroscopic cavitation bubble, are disregarded and the present work is focused on the thermodynamic effect. To account for thermal effects, it is necessary to consider in Rayleigh equation the vapor pressure corresponding to the actual temperature $T_c$ inside the bubble, which is different from the liquid temperature $T_\infty$ far from the bubble.

Using the pressure coefficient previously calculated

$$C_p = \frac{p - p_u}{(1/2)\rho_\ell W_u^2} \quad (12)$$

and the usual cavitation number defined on the basis of the vapor pressure at the liquid temperature at infinity $p_v(T_\infty)$

$$\sigma_v = \frac{p_u - p_v(T_\infty)}{(1/2)\rho_\ell W_u^2} \quad (13)$$

Rayleigh equation (11) becomes

$$R \frac{d^2 R}{dt^2} + \frac{3}{2} \left( \frac{dR}{dt} \right)^2 = \frac{p_v(T_c) - p_v(T_\infty)}{\rho_\ell} - W_u^2 \frac{C_p}{2} + \sigma_v \quad (14)$$

The first term on the right hand side accounts for thermal effects.

The bubble is supposed to travel on the blade with the local fluid velocity $W_u \sqrt{1-C_p}$. Time derivatives in the original Rayleigh equation are then transposed into space derivatives using the following equation:
$$ \frac{d}{dt} = W_u \sqrt{1 - C_p} \frac{d}{dx} \qquad (15) $$

In addition, lengths are made nondimensional using the foil chord length c. Then, Eq. (14) takes the following nondimensional form:

$$ \begin{aligned} & [1 - C_p] \left[ \bar{R} \frac{d^2 \bar{R}}{d\bar{x}^2} + \frac{3}{2} \left( \frac{d \bar{R}}{d\bar{x}} \right)^2 \right] - \frac{1}{2} \frac{d C_p}{d \bar{x}} \bar{R} \frac{d \bar{R}}{d \bar{x}} + \frac{1}{\rho_{\ell} W_u^2} \frac{d p_v}{dT} (T_{\infty} - T_c) \\ & = - \frac{C_p + \sigma_v}{2}, \end{aligned} \qquad (16) $$

where bars denote nondimensional variables.

**3.2 Convective Approach.** The problem is now to estimate the temperature inside the bubble $T_c$. Two different approaches are considered. In the present section, heat transfer at the bubble wall is supposed to be of convective nature. In other words, the heat flux at the interface is supposed to be proportional to the temperature difference $T_{\infty}-T_c$:

$$ \varphi = h(T_{\infty} - T_c) \qquad (17) $$

where h is the convection heat transfer coefficient. The validity of the convective approach is open to discussion especially when considering that the bubble is traveling at the same velocity as the liquid. From this viewpoint, the convective approach considered in Sec. 3.3 might appear more appropriate. Nevertheless, the statement that the heat flux at the bubble wall is proportional to the temperature difference appears physically reasonable. In the absence of any validated correlation applicable to cavitating flows, the heat transfer coefficient is considered here as a tuning parameter.

Under this assumption, the heat balance for the bubble is written as

$$ \frac{d}{dt} \left[ \frac{4}{3} \pi R^3 \rho_v \right] L = 4 \pi R^2 h (T_{\infty} - T_c) \qquad (18) $$

This equation allows us to compute the temperature inside the bubble:

$$ T_{\infty} - T_c = \frac{\rho_v L}{h} \frac{dR}{dt} \qquad (19) $$

Equation (19) shows that the thermal effect $\Delta T=T_{\infty}-T_c$ is proportional to $dR/dt$. In the growing phase, the bubble temperature is then lower than the liquid temperature and is higher during collapse. Using form (19) of the thermal effect, Rayleigh equation (16) becomes

$$ [1 - C_p] \left[ \bar{R} \ddot{R} + \frac{3}{2} \dot{R}^2 \right] - \frac{1}{2} \dot{C}_p \bar{R} \dot{R} + \frac{1}{\text{Nu}} \frac{\Delta p_v^*}{\rho_{\ell} \alpha_{\ell}} \tau \sqrt{1 - C_p \bar{R}} = - \frac{C_p + \sigma_v}{2} \qquad (20) $$

where Nu is the Nusselt number defined by

$$ \text{Nu} = \frac{hc}{\lambda_{\ell}} \qquad (21) $$

and $\tau=c/W_u$ is the transit time.

According to this approach, the thermodynamic effect introduces an additional term in Rayleigh equation, which is proportional to the bubble growth rate $\bar{R}$. The magnitude of thermal effects depends upon the factor:

$$ \frac{1}{\text{Nu}} \frac{\Delta p_v^*}{\rho_{\ell} \alpha_{\ell}} \tau \qquad (22) $$

It is the ratio of two characteristic times, the transit time $\tau$ and the following time:

$$ \tau_T = \text{Nu} \frac{\rho_{\ell} \alpha_{\ell}}{\Delta p_v^*} \qquad (23) $$

which can be considered as a thermal time since it essentially contains the information on heat transfer at the interface.

If the thermal time is much larger than the transit time, thermal effects are negligible and Eq. (20) reduces to the usual Rayleigh equation. Conversely, if the thermal time is much smaller than the transit time, thermal effects are predominant. A suitable criterion for estimating the magnitude of thermal effects is then based on the comparison of these two characteristic times. Besides the Nusselt number, the thermal time depends on $\Delta p_v^*/\rho_{\ell}\alpha_{\ell}$, which is very similar to Brennen's parameter (4). Equation (20) will be solved to compute the bubble evolution with thermal effects. As already mentioned, the convection heat transfer coefficient or the corresponding nondimensional Nusselt number is a free parameter of the model. We can also imagine correlating the Nusselt number with the Reynolds (Re=cW_uρ_ℓ/μ_ℓ) and Prandtl (Pr=μ_ℓ/α_ℓρ_ℓ) numbers as in traditional convection heat transfer theory, but this has not yet been done for cavitating flows.

Let us observe that Eq. (19) takes the following nondimensional form:

$$ B = \frac{\text{Pe} \sqrt{1 - C_p}}{\text{Nu} \bar{R}} \qquad (24) $$

where Pe is the Peclet number (see Eq. (27)).

**3.3 Conductive Approach.** The second approach tested here is the conductive approach. It requires to solve the heat diffusion equation in the liquid surrounding the bubble. An explicit expression of the bubble wall temperature has been proposed by Plesset and Zwick [23–25] (see also Brennen [3]) under the assumption that the thickness of the thermal boundary layer is much smaller than the bubble radius. The validity of this assumption is evaluated at the end of the present section. If so, bubble wall temperature $T_c$ can be calculated at each instant $t$ by the following expression:

$$ B = \frac{T_{\infty} - T_c}{\Delta T^*} = \frac{1}{\sqrt{\pi \varepsilon \alpha_{\ell}}} \int_{u=0}^{u=t} \frac{R^2(u) (dR/dt)(u) du}{\sqrt{\int_{v=u}^{v=t} R^4(v) dv}} \qquad (25) $$

This equation is the solution of the unsteady heat diffusion problem with a moving spherical boundary. It includes the enhancement coefficient $\varepsilon$ discussed in Sec. 1. It is similar to some extent to Eq. (19), although more complicated. Both approaches will be compared in detail in Sec. 4.1.

Using the nondimensional procedure given in Sec. 3.1 and switching to a space dependency rather than a time dependency, Eq. (25) becomes

$$ B = \frac{T_{\infty} - T_c}{\Delta T^*} = \sqrt{\frac{\text{Pe}}{\varepsilon}} J(\bar{x}) \qquad (26) $$

where Pe is the Peclet number defined by

$$ \text{Pe} = \text{Re} \text{Pr} = \frac{W_u c}{\alpha_{\ell}} \qquad (27) $$

and J is the integral

$$ J(\bar{x}) = \frac{1}{\sqrt{\pi}} \int_{u=0}^{u=\bar{x}} \frac{\bar{R}^2(u) d\bar{R}/d\bar{x}(u) du}{\sqrt{\int_{v=u}^{v=\bar{x}} \bar{R}^4(v)/\sqrt{1-C_p} dv}} \qquad (28) $$

By introducing expression (26) of the interface temperature into Rayleigh equation (16), we get the final following equation:
$$[1 - C_p] \left[ \frac{\ddot{R}}{R} + \frac{3}{2} \dot{R}^2 \right] - \frac{1}{2} \dot{C}_p \frac{\ddot{R}}{R} + \frac{\Delta p_v^*}{\rho_v W_u^2} \sqrt{\frac{\text{Pe}}{\varepsilon}} J = - \frac{C_p + \sigma_v}{2} \qquad (29)$$

When using Brennen's parameter $\Sigma$ defined in Eq. (4), previous equation takes the equivalent form:

$$[1 - C_p] \left[ \frac{\ddot{R}}{R} + \frac{3}{2} \dot{R}^2 \right] - \frac{1}{2} \dot{C}_p \frac{\ddot{R}}{R} + \Sigma \sqrt{\frac{c}{W_u^3}} \frac{1}{\sqrt{\varepsilon}} J = - \frac{C_p + \sigma_v}{2} \qquad (30)$$

Bubble evolution including thermal effects is obtained by solving Eq. (30) together with Eq. (28). The numerical procedure is presented in Sec. 3.4.

As mentioned previously, the consistency of the approach depends upon the relative order of magnitude of the thermal boundary layer thickness compared to bubble size, which is typically of the order of the chord length. The order of magnitude of the boundary layer thickness is

$$\delta \approx \sqrt{\varepsilon \alpha_{\ell} \tau} \qquad (31)$$

so that

$$\frac{\delta}{c} \approx \frac{\sqrt{\varepsilon}}{\sqrt{\text{Pe}}} = \frac{\sqrt{\varepsilon}}{\sqrt{\text{Re Pr}}} \qquad (32)$$

A typical value of the Reynolds number for rocket turbopump inducers is $10^6$. For Refrigerant 114 considered later in this work, we have $Pr \approx 5$. A typical value for the amplifying coefficient $\varepsilon$ is $10^4$ (see Sec. 4). Hence, we have $\delta/c \approx 0.04$ so that the hypothesis of a thin thermal boundary layer compared to the characteristic length scale of the cavities appears appropriate. The larger the $\varepsilon$, the larger the effective thermal diffusivity and the less accurate this hypothesis.

In a few cases and particularly when thermal effects become important, the conductive model presented above may lead to an unrealistic shape of the cavity characterized by oscillations. This effect is generally observed around cavity closure when the bubble starts to collapse, i.e., when the local pressure exceeds the vapor pressure. It is due to the fact that the temperature increases because of condensation as well as the corresponding vapor pressure. If thermal effects are large enough, the vapor pressure may exceed the local pressure. The bubble will then grow again so that the cavity interface oscillates.

To avoid this effect, which is not observed experimentally, it is chosen to set the temperature difference $T_c - T_\infty$ to zero as soon as it becomes positive. This is equivalent to assuming that the bubble temperature cannot exceed the liquid temperature. In other words, it is assumed that heat diffusion is considerably enhanced during the collapse phase so that the temperature becomes uniform almost instantaneously. The instability of the spherical shape together with the splitting of the vapor phase into small scale structures during collapse can justify this assumption. The entire growth phase is not affected by this effect. A similar assumption is made by Watanabe et al. [7], who neglected the heat flux due to condensation near the trailing edge of the cavity, assuming that condensation occurs far downstream of the cascade because of the shedding of vapor structures by the cavity.

For coherence and to make relevant comparisons between the conductive and the convective model, the same assumption of a uniform temperature in the condensation zone is also made for the convective model. The consequence on the temperature distribution in the cavity is clearly visible on Fig. 5(b), for instance.

**3.4 Numerical Aspects.** Equation (20) is a nonlinear differential equation of the second order, whereas Eq. (30) is an integrodifferential equation because of integral (28). To solve these equations, a combination of two Runge-Kutta methods of orders 4 and 5, respectively, is used. The first provides the solution, and the second an estimate of the error. This combination, due to Fehlberg

<figure><img src="image_7.png" /><figcaption>Fig. 5 Comparison between the convective (Nu=7.31×10⁶) and the conductive (ε=5×10⁴) models for the same cavity length. (a) cavity shape. (b) Temperature distribution (R114 40°C-σv=0).</figcaption></figure>

and improved by Kash and Carp [26], allows us to control the error during iterations. Indeed, for a specified tolerance, the order of magnitude of the error is calculated and the integration step is modified so that the precision criterion is satisfied: the step value is increased if the error is too small and decreased if the error is too large. Thus, the integration step is variable and the computing time is considerably reduced while keeping an excellent relative precision (10⁻⁶). This method is extended to the resolution of a system of two first order differential equations, using $R$ and $dR/dx$ as unknown functions. In this case, the relative error is simply the maximum value of the two relative errors.

The iterative resolution procedure makes it possible to transform the integrodifferential equation (30) into an ordinary differential equation. Indeed, at the current iteration, integral $J$ is calculated using estimates obtained at the previous iteration.

Mechanical and thermal equilibriums are assumed as initial conditions for the bubble. So, the interface velocity ($dR/dx)_0$ and acceleration ($d^2R/dx^2)_0$ are initially zero, and the temperature is equal to $T_\infty$ and $J_0=0$ (Eq. (25)). To correctly initiate the iterative procedure, a nonzero value of $R_0$ is required. It has been checked that the subsequent behavior of the bubble is insensitive to the value of its initial radius within a very large range of variation.

Under these assumptions, strictly speaking, Eq. (20) or (30) imposes that $C_{p0} + \sigma_v = 0$ at the initial time. The cavity starting point being fixed, this last condition is generally not satisfied, except for the very special case $\sigma_v = -C_{p0}$. Nevertheless, for all values of $\sigma_v$, the preceding initial conditions are preserved. Hence, the interface velocity presents a discontinuity at the initial time. The numerical results show that this discontinuity does not affect the resolution of the differential equation.

# 4 Results and Discussion

**4.1 Comparison Between the Convective and the Conductive Model.** Cavity shape and temperature distribution obtained by both approaches are compared on Fig. 5. Values of Nusselt number Nu for the convective model and of the ratio $\varepsilon$ of the eddy
diffusivity to the molecular one for the conductive model have been chosen so that the cavity length is the same. Figure 5(b) shows that temperature distributions are quite different for both cases.

For the convective model, temperature exhibits a sudden drop at the leading edge of the cavity. The temperature depression is maximum there and the temperature increases regularly along the cavity. As shown by Eq. (19), the convective model assumes that the temperature depression is proportional to the bubble growth rate. Since the latter is maximum at the cavity leading edge where pressure is minimum, the temperature depression is also maximum at this location.

For the conductive model, the temperature depression is zero at the cavity leading edge. The temperature inside the cavity first decreases and then increases so that a minimum is observed. This smooth behavior is the consequence of the solution of the heat diffusion equation in the liquid, which does not allow any jump in temperature.

Because of the large temperature drop at the cavity leading edge, the initial bubble growth rate is smaller for the convective approach. In the conductive case, it is the same as that without thermal effects (see Fig. 9), since temperature depression is initially zero. As a consequence, the cavity is somewhat thinner for the convective model.

The maximum cavity thickness is located closer to the cavity trailing edge for the convective approach (see Fig. 5(a)). Its location corresponds exactly to the point where the B factor vanishes, still because of the proportionality between B and $\dot{R}$. For the conductive model, the relationship between B and $\dot{R}$ is more complicated (see Eq. (25)), since the cavity temperature at any time depends upon the whole temperature history of the surrounding liquid. As a consequence, the conductive model predicts a shift between the location of maximum cavity thickness and the point where the B factor vanishes, as shown in Fig. 5(a).

It is difficult to know which model is more appropriate. On the one hand, the instantaneous initial jump in temperature predicted by the convective model may be considered as unrealistic from a physical viewpoint. On the other hand, experimental results (see Hord [6] and Fruman et al. [8]) tend to prove that the temperature depression is maximum at the cavity leading edge and decreases regularly downstream. However, the maximum observed experimentally is usually not as pronounced as the convective model predicts and measured temperature profiles are generally rather flat just downstream of the cavity detachment. On the whole, the conductive model may be preferred although no definite conclusion can be drawn from the present work. It will be shown in next section that the differences between both models decrease when thermal effects increase and that both models tend toward the same limit for large enough thermal effects.

## 4.2 Effect of Nu or ε on Thermal Effects.
Both models include a free parameter, which controls the amplitude of thermal effects. The tuning parameter is the Nusselt number for the convective approach and the ratio ε for the conductive one. Thermal effects are zero when Nu or ε are infinite and increase with decreasing Nu or ε.

The similar influence of Nu and ε can be understood by comparing both models. The heat flux to the interface is given by

$$q = h \cdot \Delta T \quad (33)$$

for the convective model, and by approximately

$$q \approx \lambda_t \cdot \frac{\Delta T}{\delta} \quad (34)$$

for the conductive one since $\Delta T/\delta$ is the order of magnitude of the temperature gradient.

The identification of forms (33) and (34) of the heat flux together with the estimate (31) of the boundary layer thickness leads to the following nondimensional relation:

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 6 Effect of Nu (convective model) or ε (conductive model) on (a) cavity length, (b) cavity shape, and (c) maximum B factor (R114 40°C-σv=0).</figcaption></figure>

$$Nu \approx \sqrt{\varepsilon Pe} \quad (35)$$

This equation suggests a link between the conductive and the convective models. To some extent, it can be expected that the influence of the Nusselt number in the convective model is similar to that of $\sqrt{\varepsilon Pe}$ for the conductive one. Figure 6(a) confirms that the variations of cavity length with Nu and $\sqrt{\varepsilon Pe}$ for the convective and the conductive model, respectively, are comparable. When thermal effects are increased, i.e., when either Nu or $\sqrt{\varepsilon Pe}$ are decreased, cavity length decreases in a similar way. Moreover, both models predict the same minimum cavity length for very large thermal effects.

This minimum cavity length is not zero. However, Fig. 6(b) indicates, and detailed computations confirm, that the cavity thickness as well as the cavity volume comes to zero when thermal effects are increased. The bubble growth rate $\dot{R}$ together with $\ddot{R}$
<figure><img src="image_12.png" /><figcaption>Fig. 7 Effect of increasing thermal effects on the distribution of the cavitation number σc based on the vapor pressure at the local cavity temperature p<sub>v</sub>(T<sub>c</sub>). Comparison with the pressure coefficient distribution (R114 40°C-σ<sub>v</sub>=0).</figcaption></figure>

also vanishes. Then, the inertial term on the left hand side of Rayleigh equation (11) tends to zero. As a consequence, the right hand side term should also vanish. This implies that, whatever the thermal model may be the vapor pressure p<sub>v</sub>(T<sub>c</sub>) based on the actual cavity temperature T<sub>c</sub> tends toward the local pressure at any point along the cavity. This is confirmed on Fig. 7, which clearly shows that the opposite of the cavitation number -σ<sub>c</sub> based on the actual cavity pressure approaches everywhere the pressure coefficient C<sub>p</sub> when thermal effects are increased. In the limit case corresponding to the suppression of cavitation by thermal effects, the temperature distribution inside the cavity is determined only by the pressure distribution on the blade and the vapor pressure curve of the fluid. In particular, the maximum value of the B factor is independent of the model, as shown by Fig. 6(c), and determined by the minimum pressure p<sub>min</sub> on the blade. The minimum value of the temperature inside the cavity T<sub>c min</sub> is determined by p<sub>v</sub>(T<sub>c min</sub>)=p<sub>min</sub>.

The existence of a minimum value of the cavity length when thermal effects become predominant can be qualitatively understood on the basis of the convective model. The limit shape of the cavity is given by the solution of Rayleigh equation (20) in which the inertial terms are neglected as previously mentioned. This simplified equation writes

$$ \dot{\bar{R}} \approx - \frac{\tau_T C_p + \sigma_v}{\tau 2 \sqrt{1 - C_p}} \quad (36) $$

Since

<figure><img src="image_13.png" /><figcaption>Fig. 8 Cavity length versus cavitation number. Comparison between computation and experiments (Ref. [22]) for cold water and R114 at two different temperatures (20 and 40°C). Cavity length is made nondimensional using blade spacing H (convective model, Nu=2.085×10<sup>6</sup>; conductive model, ε=5000).</figcaption></figure>

$$ \int_{0}^{\ell} \dot{R} dx = 0 \quad (37) $$

the minimum cavity length ℓ is determined by the following equation:

$$ \int_{0}^{\ell} \frac{C_p + \sigma_v}{\sqrt{1 - C_p}} dx = 0 \quad (38) $$

It depends only upon the pressure distribution and the cavitation number. In practice, the situation is more complicated since the B factor is forced to zero in the condensation region, where Eq. (36) is then no longer valid. The situation is even more complex in the conductive case for which equations are much more complicated. Anyway, the previous simplified argument makes the idea of a minimum cavity length quite understandable when thermal effects are increased even though an exact equation for the quantitative prediction of this minimum cavity length could not be derived.

### 4.3 Comparison With Experiments.
Figure 8 presents the evolution of cavity length with the cavitation number for water and Refrigerant 114. The cavitation number plotted here is the inducer cavitation number σ<sub>inducer</sub>. It is deduced from the cavitation number of the 2D cascade σ defined in Eq. (13) by the following relation:

$$ \frac{\sigma_{\text{inducer}}}{\sigma} \approx \frac{1}{2} \left[ \frac{r}{R} \right]^2 \quad (39) $$

where r is the local radius between hub and casing where the 2D cascade is considered and R is the inducer peripheral radius used as reference length scale for the computation of the inducer cavitation number.

In the case of water, i.e., without thermal effects, the predicted curve ℓ(σ) is in reasonable agreement with experiments. Because
<figure><img src="image_14.png" /><figcaption>**Fig. 9** Influence of thermodynamic effect on the shape of the cavity. Results for water and R114 at two different temperatures (20 and 40°C) are compared. For the three cases, the cavity length is kept constant $\lambda H=0.5$ (conductive model, $\varepsilon=5000$).</figcaption></figure>

of several shortcuts in the modeling mentioned in Sec. 1, together with the oversimplification of the real shape of the inducer blades by flat plates, discrepancies are unavoidable. However, they are not prejudicial to the present discussion, which is focused on thermodynamic effect.

The influence of thermal effects computed from both models on cavity length $ℓ(σ)$ is shown in Fig. 8 for R114 at two different temperatures (20 and 40°C). The computation was conducted with the value $\varepsilon=5000$ for the conductive model. As for the convective one, the Nusselt number was adjusted to give the same cavity length as the conductive model in the particular case $\sigma_{inducer}=0$ for R114 at 20°C. This value of the Nusselt number $(2.085 \times 10^6)$ was then kept constant for all computations including at 40°C.

The comparison of experimental and computational results shows that, with these values of $\varepsilon$ and Nu, both models lead to a reduction in cavity length due to the thermal effects of the same order of magnitude as that observed experimentally. Both models can then be considered as acceptable for the modeling of thermal effects although further computations are required for an advanced validation. In particular, the choice of the Nusselt number or the $\varepsilon$ parameter remains an open question. The use of a specific correlation between Nusselt number and Reynolds and Prandtl numbers as for classical convection would make the method predictive and give a basis for transposition of experimental results to other fluids and/or other operating conditions.

Figure 9 presents a comparison of cavity shapes obtained with the conductive model, in water (without thermal effects) and in Refrigerant 114 for two different temperatures. For each computation, the $\sigma$ value was adjusted in order to keep the cavity length approximately constant. As already observed, the cavity becomes thinner as thermal effects increase. This is due to a decrease in the bubble growth rate (and also in the collapse rate). Similarly, Ahuja and Hosangadi [16] noticed differences in the cavity content computed from their two-phase model. They report that the cavity becomes more porous and frothy as thermal effects increase as it can reasonably be expected assuming a smaller growth rate due to thermal effects.

From Fig. 9, it can be concluded that there is no similarity in cavity shape when the amplitude of thermal effects changes. Strictly speaking, cavities with and without thermodynamic effect are not similar as well as cavities for the same fluid but at two different temperatures. If thermal effects are not negligible, exact scaling requires the conservation of an additional parameter, characterizing thermal effects (as $\Sigma\sqrt{c/V^3}$ introduced in Sec. 1) in addition to the usual $\sigma$ scaling law. Results can be transposed rigorously only if such an additional scaling law is satisfied. The same conclusion holds for the temperature distribution inside the cavity which cannot be transposed from water to R114 neither between two different temperatures in R114.

The computed maximum temperature drop inside the cavity is plotted as a function of cavity length on Fig. 10 and compared to

<figure><img src="image_15.png" /><figcaption>**Fig. 10** B factor versus cavity length for $\varepsilon=5000$. Comparison between computation by the conductive model and experiments (Ref. [22]) for R114 at two different temperatures (20 and 40°C). Cavity length is made nondimensional using blade spacing H. The plotted value of B corresponds to maximum temperature depression.</figcaption></figure>

the experimental results. The experimental procedure for computing the B factor consisted in determining the equivalent vapor pressure drop (and so the corresponding temperature drop), which would lead to the observed reduction in cavity length in comparison with reference tests in cold water [22]. In the absence of any direct measurement of the temperature distribution inside the cavities, this experimental procedure based on visualizations is the only way to get a practical estimate of the temperature depression in the rotating machinery and make a possible comparison with computational results.

The predicted temperature depression increases with cavity length, as shown by experiments. The conductive approach predicts a rapid increase in the B factor for short cavities, whereas the experimental trend is an almost linear increase of the B factor with cavity length over the whole range of variation of cavity length investigated experimentally. As for the influence of the fluid temperature on thermal effects, the present approach appears to be in reasonable agreement with experiments, since it predicts a negligible influence of the fluid temperature on the nondimensional temperature depression B at least for small enough cavities as observed experimentally.

The conductive model is suitable for a comparison with experimental data in terms of B factor since the temperature distribution in the cavity is rather flat (see Fig. 5). The minimum temperature can then be considered as representative of a kind of mean cavity temperature as determined experimentally from visualizations of cavitation extent. This is not the case for the convective approach because of the sudden jump in temperature at the cavity leading edge (see Fig. 5). The temperature is far from being uniform and the comparison with experiments in terms of B factor is difficult. The mean cavity temperature, for example, proved to be inadequate for the comparison. Nevertheless, it has to be noticed that the experimental evolution of the B factor with cavity length shown in Fig. 10 was deduced from the $ℓ(σ)$ curves presented in Fig. 8 and then contains basically the same information. Hence, in the absence of direct measurements of temperature inside the cavities, it seems better to base the confrontation of computational and experimental results on the curves $ℓ(σ)$ rather than $B(ℓ)$.

## 5 Concluding Remarks

The present work is devoted to an analysis of thermal effects in a cavitating inducer. The rotating machinery is modeled by a 2D cascade of flat plates. The potential flow is computed using a
nonlinear approach based on the boundary element method. The
development of cavitation is predicted by solving Rayleigh equa-
tion using the pressure distribution deduced from the fully wetted
flow computation. The interface of the sheet cavity is approxi-
mated by the envelope of a hemispherical bubble traveling on the
suction side of the blades. Thermal effects in the Rayleigh equa-
tion are taken into account by considering the vapor pressure at
the actual bubble temperature T<sub>c</sub>, which is different from the liq-
uid temperature T<sub>∞</sub> far from the bubble.

Two different approaches are evaluated to model heat transfer
at the bubble wall and compute the cavity temperature T<sub>c</sub>. In the
convective approach, the heat flux is supposed to be proportional
to the temperature difference ΔT=T<sub>∞</sub>−T<sub>c</sub>. The amplitude of ther-
mal effects is controlled by the magnitude of the convective heat
transfer coefficient or the corresponding nondimensional Nusselt
number. The second approach is based upon the resolution of the
heat diffusion equation in the liquid surrounding the bubble. For
this conductive-type model, thermal effects are controlled by the
ratio ε of the eddy thermal diffusivity to the liquid molecular
thermal diffusivity.

The main conclusions of the present work are the following.

(1) The convective model predicts a sudden drop in tempera-
ture at the cavity leading edge followed by a gradual in-
crease up to T∞. For the conductive model, the temperature
decreases smoothly from T∞ at cavity detachment, and ex-
hibits a minimum before gradually increasing to T∞. Such
differences in temperature distributions induce differences
in cavity shapes. The convective approach predicts a thin-
ner cavity and the location of maximum thickness is also
different.

(2) Thermal effects are increased by decreasing either the Nus-
selt number for the convective model or the ε heat transfer
enhancement factor for the conductive one. Cavity length
decreases and reaches a minimum value, which is the same
for both models but which is not zero. The cavity thickness
also decreases but tends toward zero. Both models tend
toward the same solution when thermal effects are pre-
dominant. This solution corresponds to the suppression of
cavitation by thermal effects, which is achieved by a local
vapor pressure p<sub>v</sub>(T<sub>c</sub>) equal everywhere to the local pres-
sure on the blade.

(3) The evolution of cavity length with cavitation number is compared to experimental data obtained for an inducer [22] in the case of cold water (i.e., without thermal effect) and R114, which exhibits a significant thermal delay at 20 and 40°C. On the whole, both models lead to an acceptable reduction in cavity length as compared to experiments when the Nusselt number or the eddy thermal diffusivity is kept constant.

Such a simple model based on a 2D potential flow computation
coupled to the resolution of the Rayleigh equation presents clear
limitations due to simplifying assumptions. However, this work
shows that they are easy to run and can supply a support for the
interpretation of experimental data relative to thermal delays, in
particular, in inducers. They can also help in the prediction of
basic trends following changes in fluid, operating conditions and
geometry.

Acknowledgment

This research was supported by CNES/SNECMA (Contract No.
2005-015-I). The authors are particularly grateful to H. Horiguchi
and Y. Tsujimoto for their help in the validation of the potential
flow computation and to B. Chnut for his contribution to this
work.

Nomenclature

$B = B$ factor of Stepanoff (Eq. (2))

$c$ = chord length of the blade

$C_p$ = pressure coefficient

$c_{p\ell}$ = liquid heat capacity

$e$ = blade thickness

$h$ = convection heat transfer coefficient

$\ell$ = cavity length

$L$ = latent heat of vaporization

Nu = Nusselt number based on chord length (Eq. (21))

$p$ = pressure

$p_u$ = reference upstream pressure

$p_v$ = vapor pressure

Pe = Peclet number (Eq. (27))

Pr = Prandtl number

$R$ = bubble radius

$\bar{R}$ = nondimensional bubble radius $R/c$

Re = Reynolds number based on chord length

$t$ = time

$T$ = temperature

$T_c$ = local temperature in the cavity

$T_\infty$ = liquid temperature at infinity

$V$ = flow velocity

$V_u$ = upstream flow velocity

$W$ = relative flow velocity

$W_u$ = upstream relative flow velocity

$x$ = curvilinear abscissa on the blade

$\bar{x}$ = nondimensional distance $x/c$

$\alpha$ = Kato's thermodynamic parameter (Eq. (5))

$\alpha_a$ = angle of attack

$\alpha_\ell$ = thermal diffusivity of the liquid

$\alpha_t$ = turbulent thermal diffusivity

$\alpha_v$ = void fraction

$\delta$ = thermal boundary layer thickness

$\Delta p_v^*$ = characteristic vapor pressure drop (Eq. (3))

$\Delta T$ = temperature drop in the cavitating flow $T_\infty - T_c$

$\Delta T^*$ = characteristic temperature drop (Eq. (1))

$\epsilon$ = heat transfer enhancement coefficient $\alpha_t / \alpha_\ell$ or $\lambda_t / \lambda_\ell$

$\lambda_\ell$ = thermal conductivity of the liquid

$\lambda_t$ = turbulent thermal conductivity

$\mu_\ell$ = water molecular viscosity

$\rho_v$ = vapor density

$\rho_\ell$ = liquid density

$\Sigma$ = Brennen's thermodynamic parameter (Eq. (4))

$\sigma_c$ = local cavitation number based on $p_v(T_c)$

$\sigma_v$ = cavitation number (Eq. (10))

$\tau$ = transit time $c/W$

References

[1] Stahl, H. A., Stepanoff, A. J., and Phillipsburg, N. J., 1956, "Thermodynamic Aspects of Cavitation in Centrifugal Pumps," ASME J. Basic Eng., November 1956, pp. 1691-1693.

[2] Stepanoff, A. J., 1964, "Cavitation Properties of Liquids," J. Eng. Power, pp. 195-200.

[3] Brennen, C. E., 1995, *Cavitation and Bubble Dynamics*, Oxford University Press, New York.

[4] Brennen, C. E., 1994, *Hydrodynamics of Pumps*, Oxford University Press, New York.

[5] Kato, H., 1984, "Thermodynamic Effect on Incipient and Developed Sheet Cavitation," International Symposium on Cavitation Inception, New Orleans, Dec. 9-14, Vol. 16, pp. 127-136.

[6] Hord, J., "Cavitation in Liquid Cryogens," NASA Report No. CR-2054; 1972, NASA Report No. CR-2156; 1973, NASA Report No. CR-2242; 1974, NASA Report No. CR-2448.

[7] Watanabe, S., Hidaka, T., Horiguchi, H., Furukawa, A., and Tsujimoto, Y., 2005, "Steady Analysis of Thermodynamic Effect of Partial Cavitation Using Singularity Method," *2005 ASME Fluids Engineering Division Summer Meeting and Exhibition*, Houston, TX, Jun. 19-23.

[8] Fruman, D. H., Benmansour, I., and Sery, R., 1991, "Estimation of the Thermal Effects on Cavitation of Cryogenic Liquids," 1991, *Cavitation and Multiphase Flow Forum*, ASME FED, 109, pp. 93-96.
[9] Fruman, D. H., Reboud, J. L., and Stutz, B., 1999, "Estimation of Thermal Effects in Eavitation of Thermosensible Liquids," Int. J. Heat Mass Transfer, 42, pp. 3195-3204.

[10] Billet, M. L., 1970, "Thermodynamic Effects on Developed Cavitation in Water and Freon 113," MS thesis, The Pennsylvania State University.

[11] Billet, M. L., Holl, J. W., and Weir, D. S., 1981, "Correlations of Thermodynamic Effects for Developed Cavitation," J. Fluids Eng., 103, pp. 534-542.

[12] Holl, J. W., Billet, M. L., and Weir, D. S., 1975, "Thermodynamic Effects on Developed Cavitation, J. Fluids Eng., 97, pp. 507-514.

[13] Rapposelli, E., and d'Agostino, L., 2003, "A Barotropic Cavitation Model with Thermodynamic Effects, Fifth International Symposium on Cavitation, Osaka, Japan, Nov. 1-4.

[14] Rolland, J., Boitel, G., Barre, S., Goncalves, E., and Fortes, Patella R., 2006, "Experiments and Modelling of Cavitating Flows in Venturi-Part I: Stable Cavitation," Sixth International Symposium on Cavitation, Wageningen, The Netherlands, Sept.

[15] Tani, N., and Nagashima, T., "Numerical Analysis of Cryogenic Cavitating Flow on Hydrofoil: Comparison Between Water and Cryogenic Fluids," Fourth International Contevence on Launcher Technology, Liege, Belgium, Dec. 3-6.

[16] Ahuja, V., and Hosangadi, A., 2006, "A Numerical Study of Cavitation in Cryogenic Fluids. Part I: Mean Flow Parametric Studies," Sixth International Symposium on Cavitation, Wageningen, The Netherlands, Sept.

[17] Hosangadi, A., and Ahuja, V., 2006, "A Numerical Study of Cavitation in Cryogenic Fluids. Part II: New unsteady Model for Dense Cloud Cavitation,"

Sixth International Symposium on Cavitation, Wageningen, The Netherlands, Sept.

[18] Hosangadi, A., and Ahuja, V., 2005, "Numerical Study of Cavitation in Cryogenic Fluids," J. Fluids Eng., 127, pp. 267-281.

[19] Kubota, A., Kato, H., and Yamaguchi, H., 1992, "A New Modeling of Cavitating Flows: A Numerical Study of Unsteady Cavitation on a Hydrofoil Section," J. Fluid Mech., 240, pp. 59-96.

[20] Singhal, A. K., Athavale, M. M., Li, H., and Jiang, Y., 2002, "Mathematical Basis and Validation of the Full Cavitation Model," J. Fluids Eng., 124, pp. 617-624.

[21] Franc, J. P., and Michel, J. M, 2004, Fundamentals of Cavitation, Kluwer Academic, Dordrech; Franc, J. P., Janson, E., Morel, P., Rebattet, C., and Riondet, M., 2001, "Visualizations of Leading Edge Cavitation in an Inducer at Different Temperatures," Fourth International Symposium on Cavitation, Pasadena, CA, Jun. 20-23.

[22] Franc, J. P., Rebattet, C., and Coulon, A., 2004, "An Experimental Investigation of Thermal Effects In a Cavitating Inducer," J. Fluids Eng., 126, pp. 716-723.

[23] Plesset, M. S., and Zwick, S. A., 1952, "A Nonsteady Heat Diffusion Problem with Spherical Symmetry," J. Appl. Phys., 23(1), pp. 95-98.

[24] Plesset, M. S., and Zwick, S. A., 1954, "The Growth of Vapour Bubbles in Superheated Liquids," J. Appl. Phys., 25(4), pp. 493-500.

[25] Zwick, S. A., and Plesset, M. S., 1955, "On the Dynamics of Small Vapour Bubbles in Liquids," J. Math. Phys. (Cambridge, Mass.), 33, pp. 308-330.

[26] Fortin, A, 2002, Analyse numérique, Internationales Polytechniques, Montreal.