Numerical simulation of cavitating flows under uncertainty

Cavitation is characterized by vapor bubbles creation in the liquid phase as a consequence of a pressure drop. This phenomenon can be reproduced by means of several two-phase models. An equation of state is commonly used in order to define the thermophysical properties of the two fluids and to close the model. The aim of this work is to study how the uncertain parameters of the equation of state (EOS) can influence the prediction of the cavitation structures. These uncertainties are propagated through a two-phase numerical solver for evaluating the impact on the predictive character of the numerical solution. The variability of the mixture velocity and the mixture pressure are analyzed.


Introduction
The cavitation is a phenomenon affecting several mechanical devices. Classical examples are the propellers, the hypervelocity projectiles and all devices in which a section decrease can generate a very fast liquid acceleration, as in the pump or in the rocket engines. Since the nucleation of a new phase appears, cavitation can yields dramatic consequences, such as material erosion [21], noise [22] and performance decrease [23].
Several numerical models have been proposed to simulate this phenomenon. Principally two classes can be identified: (1) methods considering interfaces as true discontinuities and (2) methods allowing a numerical diffusion at the interface. Methods belonging to the second class treat the interface like an artificial transition region where the thermodynamic conditions are unknown and they are known as Interface diffusive Models [5,4,3].
In this work, a five equation model, with heat and mass transfer terms, is used for reproducing cavitating flows, supposing that the two phases have the same velocity and pressure. The numerical approach proposed in [6] has been followed here. It is based on a splitting procedure based on three steps for assuring a good accuracy and a reduced computational cost. The step associated to definition of the heat and mass transfer terms is strongly connected to the EOS chosen in order to close the system. This means that the mass transfer prediction is sensitive to the choice of the EOS parameters [2,7,6].
Despite the use of stochastic methods applied to the numerical simulation in fluid mechanics being more and more diffused, only few studies exist concerning the application of uncertainty quantification tools to cavitating flows. Li et al. [9] proposed a Markov stochastic model to reproduce the random behaviour of cavitation bubble(s) near compliant walls. Fariborza et al. [8] proposed an empirical model for the time-discrete stochastic nucleation of intergranular creep cavities. They assumed nucleation to occur randomly in time, with the temporal behavior being governed by an inhomogeneous Poisson process. Giannadakis et al. [10] described the bubble breakup in lagrangian models using a stochastic Monte-Carlo approximation. This study was oriented on the particular topic of cavitation in the Diesel nozzle holes. Wilczynski [11] and Goel et al. [12] performed an uncertainties-based study on some hydrodynamic cavitation model parameters. In particular, Wilczynski [11] applied a stochastic model to capture the interaction of turbulent pressure field on cavitation nuclei population. Moreover, Goel et al. [12] performed a sensitivity analysis on several empirical parameters used typically in two-phase models. This study was performed using a finite differences method. In this case, input data uncertainty characterization is not required for the sensitivity analysis, that can be performed basing only on the mathematical form of the model. A recent work [13] presents a sensitivity study of the cavitation and turbulence models, where different combinations of empirical coefficients are considered. Recently, Rodio & Congedo [14,26,25] have proposed a study about the impact of various sources of uncertainty (on the cavitation model and on the inlet conditions) on the prediction of cavitating flows by coupling a non-intrusive Polynomial Chaos stochastic method with a cavitating CFD solver. The aim of this paper is to evaluate the influence of experimental and EOS uncertainties on the cavitation model prediction. Note that a Stiffened Gas EOS is used here.
First, the uncertainties and their variation are estimated. Then a forward problem is applied, i.e. the uncertainties are propagated through a two-phase numerical solver for evaluating the impact on the predictive character of the numerical solution. In particular, the variability of some quantities of interest, such as the cavitation length, mixture pressure, is computed, thus permitting to evaluate the robustness of the physical model.
The paper is organized as follows. Section 2 presents the five-equation model with heat and mass transfer. Section 3 illustrates the main ingredients of the splitting method used for solving the cavitation model. In Section 4, the Stiffened Gas equation of state is be described. In section 5, some details concerning the Uncertainty Quantification method are provided. Finally, section 6 illustrate numerical results obtained on the most used numerical configuration in literature.

The cavitation model
Supposing a mechanical relaxation between the two phases, i.e. supposing that the two phases have the same pressure and velocity, a five equations model can be used for reproducing a cavitation problem. It is composed by a transport equation for the vapor volume fraction, the mass conservative equations for each phase and, finally, the momentum and the energy conservative equation for the liquid/vapor mixture, as follows: where α 1 , ρ k , v, P , e k , Q andẎ are the gas volume fraction, the phase density, the mixture velocity, the mixture pressure, the specific phase internal energy , the heat term and the mass term, respectively. E = e + 1 2 v 2 is the mixture total energy and e = (α 1 ρ 1 e 1 + α 2 ρ 2 e 2 )/ρ is the mixture internal energy. The parameters K = . For closing the system, an equation of state is required, in order to define e k as a function of pressure and density for both phases, and the definition of Q andẎ . This is detailed in Section 3.

A splitting method
The system (1) is solved by applying a splitting method proposed in [6]. Assuming that the characteristic time of a mechanical relaxation, 1/µ, is much smaller than the characteristic time scales 1/θ and 1/ν of heat and mass transfer terms, respectively, the solution at each time can be found by performing the following three steps : • Step 1 : compute the numerical solution of the hyperbolic part of system (1) without heat and mass transfer term source. • Step 2 : update the solution of system (1), by solving the temporal ODEs system with the heat and mass transfer terms. This consists in applying a thermo-chemical relaxation that allows to attain a new equilibrium solution on the saturation curve. • Step 3 : the thermo-chemical relaxation does not guarantee that the positivity of the solution be preserved, so a positivity check of the solution is performed.
Since the originality of this work is focused on the uncertainty quantification analysis, we will describe, in the follow-up of this section, the main ingredients of the splitting method. We suggest, anyway, to read [6] for a more accurate description of the numerical method.

Step 1: Numerical solution of a five-equation model without heat and mass transfer
This step is necessary for solving the system (1) without heat and mass transfer terms. The problem is, thus, reduced to apply any numerical method able to solve a classical five equation model, as for example in [7,3]. In this work, we applied the discrete equation method (DEM) [16], since it allows to treat the approximation of a non-conservative term, when the velocity v and K are simultaneously discontinuous [15,16,4].

3.2.
Step 2: Numerical solution of the temporal ODEs with heat and mass transfer terms This step consists in solving the following temporal ODEs, considering the new source terms: The mass and heat transfer are two unknowns of the system. Applying a thermochemical relaxation, i.e. imposing temperature and Gibbs free energy equilibrium between the phases, a non trivial manipulation of the system 2 allows to find the final expressions for Q andẎ : where: and ξ = ( , the interface density has been defined in [2] as follows: ρ I = ( In all equations T k , S k , Cv k , c k and ρ k are the phase temperature, entropy, heat capacity at constant volume, speed of sound and density respectively.
Note that the expressions of Q andẎ have been obtained supposing the use of a Stiffened Gas as equation of state. So the terms 4 are clearly dependent on the EOS parameters. Their influence on the simulation is shown in section 4.

Step 3: Mass fraction and density positivity
The previous step allows to compute the heat and mass transfer terms, but it does not guarantee a consistent solution. So, in order to preserve the positivity of the solution, a limitation is put on the source terms, by determining the maximum admissible values, as follows: Thus, if |S max,α 1 | > |S α 1 | and |S max,Y 1 | > |S Y 1 |, the source terms is used. Otherwise, if the limit value is not respected, the gas volume fraction α 1 = 1 is imposed if S 1 > S 2 , otherwise α 1 = 0. Knowing α 1 , we can find the new pressure, P, by solving a single quadratic equation (see [6] for more details): where: Then, by solving the single quadratic equation of P ⋆⋆ , we select the physically admissible solution of the quadratic equation that maximizes the total entropy s ⋆⋆ = Y ⋆⋆ /ρ ⋆⋆ This procedure allows to reduce the computational cost compared to [2,17]: it is not necessary to integrate over a fractional hydrodynamic time step.

Thermodynamic closure
The equation of state chosen for this analysis is the Stiffened Gas (SG). This EOS has been largely used for simulating two-phase flows problems. On the contrary of more complex EOSs, the expressions for each thermodynamic variables are simple and defined as follows: where P k , ρ k and e k are the pressure, the density and the specific internal energy of the phase, respectively. The polytropic coefficient γ k is the constant ratio of specific heat capacities γ k = c p,k /c v,k , P ∞,k is a constant reference pressure,q k is the specific creation energy of the fluid at a given reference state and q ′ is a fluid specific constant. Moreover, T k and c v,k are the temperature, the specific heat at constant volume and the enthalpy, respectively. Moreover, the sound speed remains strictly positive, ensuring the hyperbolicity of the system. As explained in section 3.2, the mass term is dependent on the EOS parameters. Let focus the attention on the parameter q ′ . It is usually estimated by means the procedure described in [2]. However, for the same fluid and the same test case, several values are used in literature. We focus the attention for example on the water and its vapor parameters used for the two-phase expansion tube test (see [2]). In [2,17], q ′ = −23.0 × 10 3 , in [7] q ′ = −23.4 × 10 3 and in [6] The aim of this paper is, thus, to evaluate and quantify the influence of this parameter on the model prediction.

Forward Uncertainty Quantification method
For the uncertainty propagation, a truncated Polynomial Chaos expansion (see [20]) is computed. Using this non-intrusive uncertainty quantification method means that a single deterministic computation is replaced with a whole set of computations, each one of those being run for specific values of the uncertain conditions. Polynomial Chaos (PC) expansions are derived from the original theory of Wiener on spectral representation of stochastic processes using Gaussian random variables. Let ξ be a vector of standard independent random variables ξ i , i = 1, 2, ..., n ξ . Any well-behaved process u (i.e. a second-order process, then with a finite variance) can be expanded in a convergent (in the mean square sense) series of the form where α are multi-indices, α = (α 1 , α 2 , ..., α n ), with each component α i = 0, 1, ..., and Ψ α are multivariate polynomial functions orthogonal with respect to the probability distribution function of the vector ξ. Each Ψ α is defined by a product of orthogonal polynomials A one-to-one correspondence exists between the choice of stochastic variable ξ i and the polynomials Φ α i i (ξ i ). For instance, if ξ i is a normal/uniform variable, the corresponding Φ α i i (ξ i ) are Hermite/Legendre polynomials of degree α i . Coefficients u α (x, t) are called the PC coefficients of the random process u and are obtained by where the scalar product is defined by the expectation operator. For practical use, the PC expansions are truncated to degree N o The number of multivariate polynomials Ψ α , that is, the dimension of the expansion basis, is related to the stochastic dimension n ξ and the degree N o of polynomials ; it is given by the formula (n ξ + N o)!/(n ξ !N o!). Several approaches can be used to estimate PC coefficients. The approach used in this study is based on Gauss-Legendre quadrature nodes. From the PC expansion of the random process, it is then easy to derive its mean and variance and to estimate sensitivity information using the analysis of variance (ANOVA) decomposition [18].

Results
In this section, we present some results. The test case reproduced for the validation of the cavitation model and the uncertainties propagation is a two-phase expansion tube.
Physical parameters are summarized in Table 1. A second order scheme has been performed for all the numerical test-cases by means of a MUSCL scheme [15], coupled with a Relaxation solver [16] and a Van Leer limiter. Note also that a mesh-convergence study has been done, even if it is not reported here for brevity.

Two-phase expansion tube test
Now, let us consider a tube of unit length, filled with water at atmospheric pressure P = 10 5 , density ρ 2 = 1150 kg/m 3 and temperature T 2 = 354.728 K. A small amount of vapor, α 1 = 10 −2 , is uniformly distributed in the whole domain. By assuming the flow in thermal and pressure initial equilibrium T 2 = T 1 , the vapor density can be computed by means of the equations (7) x [m] Supposing theoretically that the liquid phase can not evaporate, because of the incremental velocity, the pressure drops quite to zero in the center of the tube, as we can observe in figure  1 (solid line). In real conditions, the pressure can drops until it attains the saturation value corresponding to the flow temperature. So, between x = 0.3 m and x = 0.7 m, the pressure reaches its saturation value of about 5 × 10 4 Pa.
The comparison with the results obtained by Pelanti and Shyue [7] and by Zein et al. [17], shows a very good agreement in terms of pressure and velocity (Fig. 1b-c), but some differences in terms of vapor volume fraction (see Fig. 1a) can be observed in particular with [17]. Note that this difference could be explained by the different value used for the equation of state's parameter q ′ .

Uncertainty assessment
In this case, three uncertainties are considered: the EOS parameter, q ′ , the initial left velocity and the vapour volume fraction. Because of the lack of knowledge concerning the parameter q ′ , it is assumed as an epistemic uncertainty. Its variation range is chosen relying on a calibration with respect to the National Institute of Standards and Technology (NIST [24]) experimental data. In particular, the saturation curves data obtained in terms of pressure against temperature and in terms of the liquid volume against temperature are considered. By varying q ′ several curves are obtained (Figures 2(a) and 2(b)). Considering the experimental uncertainty provided by NIST, q ′ is assumed to vary between -23300 and -23500. s The initial left velocity is assumed as an epistemic uncertainty, with a arbitrary variation of ±10% with respect to the initial value of -2 m/s (treated with an uniform probability density function).
The variation of the initial vapor volume fraction is established by considering the existing literature about the existence of non-condensable gas in the liquid phase. In particular, relying on experimental observations for water at room temperature [1], the initial gas volume fraction is assumed to vary between 5 · 10 −3 and 1 · 10 −2 .

Forward propagation
First, the uncertainties propagation study has permitted to compute the mean and the variance of the pressure and the velocity, obtained at the final time (t=3.2 ms) along the x coordinate. The mean profiles of the pressure is very similar to the deterministic solution, shown in figure 1b. However the numerical errors bars show a large variability of the pressure. In order to identify which is the most important uncertainty that produces this variability, an ANOVA analysis is then performed. Figure 3(b) shows the pressure variance and all the contributions to the variance related with the three uncertainties. The pressure variation is only dependent on the uncertainties on α GAS and q ′ . In particular the initial α variation produces a non-negligible pressure variation, but it does not affect the transition phase region (0.3 < x < 0.7 m), which, on the contrary, is strongly influenced by the q ′ variation.
In terms of velocity, the variation is less important respect to the pressure one (see figure  4(a)). Anyway, we remember that the velocity variation has been applied to the left part of the tube. It explains the variance behavior for x < 0.5 m. In fact, the contribution of initial velocity contribution, v IN , is more important than the others, but remains of the same order of magnitude. For x > 0.5 m, on the contrary, the v IN contribution drops quite to zero and the q ′ contribution remains the most important in the cavitation region.

Conclusions and Perspectives
This paper is focused on an uncertainty propagation study aimed to estimate the influence of experimental uncertainties on the cavitation-driven phenomena. The influence of uncertainties on boundaries condition and EOS uncertainties have been assessed on the mixture pressure and velocity. The study showed that the cavitation region is strongly influenced by the EOS parameter uncertainties. This clearly questions the quality of the prediction in this kind of configuration, considering the ranges of variation.