Simulations of unsteady cavitating turbulent flow in a Francis turbine using the RANS method and the improved mixture model of two-phase flows

Abstract

This paper reports the simulation results for the unsteady cavitating turbulent flow in a Francis turbine using the mixture model for cavity–liquid two-phase flows. The RNG k–ε turbulence model is employed in the Reynolds averaged Navier–Stokes equations in this study. In the mixture model, an improved expression for the mass transfer is employed which is based on evaporation and condensation mechanisms with considering the effects of the non-dissolved gas, the turbulence, the tension of interface at cavity and the effect of phase change rate and so on. The computing domain includes the guide vanes, the runner, and the draft tube, which is discretized with a full three-dimensional mesh system of unstructured tetrahedral shapes. The finite volume method is used to solve the governing equations of the mixture model and a full coupled method is combined into the algorithm to accelerate the solution. The computing results with the mixture model have been compared with those by the single-phase flow model as well as the experimental data. The simulation results show that the cavitating flow computation based on the improved mixture model agrees much better with experimental data than that by the single-phase flow calculation, in terms of the amplitude and dominated frequency of the pressure fluctuation. It is also observed from the present simulations that the amplitude of the pressure fluctuation at small flow rate is larger than that at large flow rate, which accords with the experimental data.

Appendices

Appendix 1

1.1 Cavitation model

1.1.1 Original cavitation mass transfer model

The model presented in the previous papers by Liu et al. [24, 25]:

The mass transfer rate S is obtained as follows:

If $p > p_{\text{v}}^{ * }$:

$$ S = {\frac{{3\alpha_{\text{v}} }}{r}}C_{1} \left( {{\frac{M}{{ 2\pi {\text{R}}}}}} \right)^{{\frac{1}{2}}} \left( {{\frac{{p_{\text{v}}^{ * } }}{\sqrt T }} - {\frac{p}{\sqrt T }}} \right). $$

(19)

If $p < p_{\text{v}}^{ * }$:

$$ S = {\frac{{3(1 - \alpha_{\text{v}} - \alpha_{u} )}}{r}}C_{2} \left( {{\frac{M}{{ 2\pi {\text{R}}}}}} \right)^{{\frac{1}{2}}} \left( {{\frac{{p_{\text{v}}^{ * } }}{\sqrt T }} - {\frac{p}{\sqrt T }}} \right) $$

(20)

with

$$ p_{\text{v}}^{ * } = p_{\text{v}} + 0.195\rho k. $$

(21)

In this model, the cavitation mass transfer is obtained by using the kinetic theory of mass transfer and is based on the mechanisms of evaporation and condensation. The cavitation mass transfer in this model is directly proportional to the difference between the static pressure p and the vapor pressure p _v. The influence of liquid surface tension λ is considered through (21) and the influence of liquid turbulence kinetic energy k is considered through equation (22). This model is the combination of both the mass flux model obtained by Cammenga [4] and the full cavitation model by Singhal et al. [35]. The cavitation mass transfer model by Singhal et al. [35] is based on the Rayleigh–Plesset bubble dynamic equation, and is directly proportional to the square root of the difference between the static pressure and the vapor pressure.

Introducing the coefficient ξ of evaporation and condensation of water, and combining (19) and (20), we obtain,

$$ S = {\frac{{3\alpha_{\text{v}} }}{r}}C_{1} {\frac{2\xi }{2 - \xi }}\left( {{\frac{M}{ 2\pi R}}} \right)^{{\frac{1}{2}}} \left( {{\frac{{p_{\text{v}}^{ * } }}{\sqrt T }} - {\frac{p}{\sqrt T }}} \right) $$

(22)

where ξ = 0.051 [6].

The vapor volume fraction α _v in (19) (or 22) and (20) can be related to the bubble number density n and radius of bubble r by $\alpha_{\text{v}} = n\frac{4}{3}\pi r^{3} .$ However, the bubble number density n is unknown, or is difficult to be determined.

1.1.2 Corrections to the cavitation mass transfer rate

For simplicity, the typical bubble size r is taken to be the same as the limiting (maximum possible) bubble size. Then, r is determined by the balance between aerodynamic drag and surface tension force. A commonly used correlation in the nuclear industry is [35],

$$ r = {\frac{0.061We\lambda }{{2\rho_{l} u_{\text{rel}}^{2} }}} $$

(23)

where We is the Webber number, u _rel is the relative speed (slipping speed) between liquid phase and cavity phase.

In bubbly flow regime, u _rel is generally fairly small, e.g., 5–10% of liquid velocity. According to Singhal et al. [35] by using various limiting arguments, e.g., r → 0 as α → 0, and to the fact that the phase change rate per unit volume should be proportional to the volume fractions of the donor phase, the cavitation transfer term can be expressed in the form:$p_{\text{v}}^{*} > p_{\text{v}}$:

$$ S = C_{1} {\frac{{3\alpha_{\text{v}} \rho_{l} \sqrt k }}{\lambda }}\left( {{\frac{2\xi }{2 - \xi }}} \right)\left( {{\frac{M}{{ 2\pi {\text{R}}}}}} \right)^{{\frac{1}{2}}} \left( {{\frac{{p_{\text{v}}^{*} }}{\sqrt T }} - {\frac{p}{\sqrt T }}} \right) $$

(24)

$p_{\text{v}}^{*} < p_{\text{v}}$:

$$ S = C_{2} {\frac{{3\left( {1 - \alpha_{u} - \alpha_{\text{v}} } \right)\rho_{l} \sqrt k }}{\lambda }}\left( {{\frac{2\xi }{2 - \xi }}} \right)\left( {{\frac{M}{ 2\pi R}}} \right)^{{\frac{1}{2}}} \left( {{\frac{{p_{\text{v}}^{*} }}{\sqrt T }} - {\frac{p}{\sqrt T }}} \right). $$

(25)

The set of (24) and (25) is the improved expression for cavitation mass transfer employed in present paper.

1.1.3 Validation of the two-phase cavitation model with mass transfer

Using the two-phase cavitation model with mass transfer used in this paper, a 3D cavity–liquid turbulent two-phase flow around a hydrofoil was calculated. As shown below, the results yielded agreements with the test data for two cases published in the literature.

The first hydrofoil is of an asymmetric leading edge (ALE25) as shown in Fig. 16. The velocity distributions predicted with present calculation model at 4 positions show that good agreement with the experiment is achieved (Fig. 17). In Fig. 17, the velocity at the inlet is 13 m/s. The pressure at the outlet is 197,200 Pa. Non-slip condition is applied on solid walls. In the present computation, the coefficients C ₁ = 0.13 (evaporation) and C ₂ = 0.01 (condensation), proposed by as Hagen et al. [16] are used.

Fig. 16

Sketch of ALE25 hydrofoil

Fig. 17

Comparison of velocity distributions on ALE25 hydrofoil between calculations and experiments. a X = 11 mm, b X = 22 mm, c X = 33 mm, d X = 44 mm

The second simulation is for the ALE15 hydrofoil [11]. The geometry of the hydrofoil and the flow conditions are similar to those shown for ALE25 above. The simulation results with the present two-phase cavitation model with mass transfer and the full cavitation model (Singhal et al. [35] from the Rayleigh–Plesset equation) are compared with the experiments (Fig. 18). The velocity distribution in y direction on the section of z = 5 mm and x = 13 mm of the hydrofoil is shown. It can be seen from Fig. 18 that the simulation result from the present two-phase cavitation model in this paper is better than that from the full cavitation model.

Fig. 18

Velocity distribution along x direction on the section of z = 5 mm, x = 13 mm of ALE15 hydrofoil

Appendix 2

2.1 Fully coupled method

The Fully coupled method is used to solve Navier–Stokes equations with primary variables of velocity components u, v, w and pressure p at all grid notes of the whole calculation domain by using one algebraic equation group at one-time step.

Figure 1 shows the control volume of a node and the interfaces. The discretized continuity equation and momentum equations in the control volume can be expressed as follows,

$$ \begin{gathered} a_{c}^{{u_{i - 1/2,j,k} }} u_{i - 1/2,j,k} + a_{c}^{{u_{i + 1/2,j,k} }} u_{i + 1/2,j,k} + a_{c}^{{v_{i,j - 1/2,k} }} v_{i,j - 1/2,k} \hfill \\ + a_{c}^{{v_{i,j + 1/2,k} }} v_{i,j + 1/2,k} + a_{c}^{{w_{i,j,k - 1/2} }} w_{i,j,k - 1/2} + a_{c}^{{w_{i,j,k + 1/2} }} w_{i,j,k + 1/2} = s_{i,j,k}^{c} \hfill \\ \end{gathered} $$

(26)

$$ a_{\text{m}}^{{u_{i - 1/2,j,k} }} u_{i - 1/2,j,k} + a_{p}^{{u_{i - 1/2,j,k} }} p_{i,j,k} = s_{\text{m}}^{{u_{i - 1/2,j,k} }} $$

(27)

$$ a_{\text{m}}^{{u_{i + 1/2,j,k} }} u_{i + 1/2,j,k} \, + a_{p}^{{u_{i + 1/2,j,k} }} p_{i,j,k} = s_{\text{m}}^{{u_{i + 1/2,j,k} }} $$

(28)

$$ a_{\text{m}}^{{v_{i,j - 1/2,k} }} v_{i,j - 1/2,k} + a_{p}^{{v_{i,j - 1/2,k} }} p_{i,j,k} = s_{\text{m}}^{{v_{i,j - 1/2,k} }} $$

(29)

$$ a_{\text{m}}^{{v_{i,j + 1/2,k} }} v_{i.j + 1/2,k} + a_{p}^{{v_{i,j + 1/2,k} }} p_{i,j,k} = s_{\text{m}}^{{v_{i,j + 1/2,k} }} $$

(30)

$$ a_{\text{m}}^{{w_{i,j,k - 1/2} }} w_{i,j,k - 1/2} + a_{p}^{{w_{i,j,k - 1/2} }} p_{i,j,k} = s_{\text{m}}^{{w_{i,j,k - 1/2} }} $$

(31)

$$ a_{\text{m}}^{{w_{i,j,k + 1/2} }} w_{i,j,k + 1/2} + a_{p}^{{w_{i,j,k + 1/2} }} p_{i,j,k} = s_{\text{m}}^{{w_{i,j,k + 1/2} }} $$

(32)

where a _i is coefficient; s _i is the source term. The group of the local linear algebraic equation may be expressed in the matrix form:

$$ \left[ {\begin{array}{*{20}c} {a_{\text{m}}^{{u_{i - 1/2,j,k} }} } & 0 & 0 & 0 & 0 & 0 & {a_{p}^{{u_{i - 1/2,j,k} }} } \\ 0 & {a_{\text{m}}^{{u_{i + 1/2,j,k} }} } & 0 & 0 & 0 & 0 & {a_{p}^{{u_{i + 1/2,j,k} }} } \\ 0 & 0 & {a_{\text{m}}^{{v_{i,j - 1/2,k} }} } & 0 & 0 & 0 & {a_{p}^{{v_{i,j - 1/2,k} }} } \\ 0 & 0 & 0 & {a_{\text{m}}^{{v_{i,j + 1/2,k} }} } & 0 & 0 & {a_{p}^{{v_{i,j + 1/2,k} }} } \\ 0 & 0 & 0 & 0 & {a_{\text{m}}^{{w_{i,j,k - 1/2} }} } & 0 & {a_{p}^{{w_{i,j,k - 1/2} }} } \\ 0 & 0 & 0 & 0 & 0 & {a_{\text{m}}^{{w_{i,j,k + 1/2} }} } & {a_{p}^{{w_{i,j,k + 1/2} }} } \\ {a_{c}^{{u_{i - 1/2,j,k} }} } & {a_{c}^{{u_{i + 1/2,j,k} }} } & {a_{c}^{{v_{i,j - 1/2,k} }} } & {a_{c}^{{v_{i,j + 1/2,k} }} } & {a_{c}^{{w_{i,j,k - 1/2} }} } & {a_{c}^{{w_{i,j,k + 1/2} }} } & 0 \\ \end{array} } \right]\left\{ {\begin{array}{*{20}c} {u_{i - 1/2,j,k} } \\ {u_{i + 1/2,j,k} } \\ {v_{i,j - 1/2,k} } \\ {v_{i,j + 1/2,k} } \\ {w_{i,j,k - 1/2} } \\ {w_{i,j,k + 1/2} } \\ {p_{i,j,k} } \\ \end{array} } \right\} = \left\{ {\begin{array}{*{20}c} {s_{m}^{{u_{i - 1/2,j,k} }} } \\ {s_{m}^{{u_{i + 1/2,j,k} }} } \\ {s_{m}^{{v_{i,j - 1/2,k} }} } \\ {s_{m}^{{v_{i,j + 1/2,k} }} } \\ {s_{m}^{{w_{i,j,k - 1/2} }} } \\ {s_{m}^{{w_{i,j,k + 1/2} }} } \\ {s^{c}_{i,j,k} } \\ \end{array} } \right\} $$

(33)

Then, all groups of the local algebraic equations at all volumes in calculated domain are assembled together. Thus, the group of the algebraic equations with 3G + (G-M) number of equations is obtained, where G is the number of total nodes at which the velocity components store; G-M is the number of total control volumes, at which the pressure values store.

In the solution of the cavitating turbulent flow, the velocity and pressure at whole domain are solved first. Then, the k and ε equations, and the bubble volume fraction equation are solved sequentially. After all the equations are solved, the iterative procedure moved to next time step.