Pietro Marco Congedo

Eric Goncalves Da Silva

Maria Giovanna Rodio

E Goncalves

About the uncertainty quantification of turbulence and cavitation models in cavitating flows simulations

Keywords: cavitation, turbulence, uncertainty quantification, Bayesian-based methods

Introduction

The collapse of vapour bubbles in high-pressure region can produce dramatic effects, such as failure, erosion and other undesirable effects. In particular, an accurate description of interactions between the vapour and liquid phases requires accurate physical models and a way to take into account the dynamics of the interface, multiscale effects, turbulence and thermodynamics.

Several numerical approaches have been proposed to reproduce cavitating flows in external and internal configurations. Principally the models can be regrouped in two major categories: interface models and two-phase models. In the first case, the liquid and the vapour phase are separated by an interface, then the systematic reconstruction of interface and the applicability to complex geometries are the most challenging issues. Concerning two-phase models, the two phases are treated as a mixture. Difficulties of these models are related to the mixture's properties estimation based on the liquid-vapour mixture ratios. Differences between the various models in the second category mostly come from the relation that defines the density field. For more details concerning the various modelling approaches, Refs. [START_REF] Agostino | Fluid Dynamics of Cavitation and Cavitating Turbopumps[END_REF][START_REF] Brennen | Cavitation and Bubble Dynamic[END_REF][START_REF] Franc | Analysis of thermal effects in a cavitating inducer using rayleigh equation[END_REF][START_REF] Utturkar | Recent progress in modeling of cryogenic cavitation for liquid rocket propulsion[END_REF] are strongly recommended.

An open question is still related to the interaction between cavitation and turbulence, in order to assess the influence of sheet cavitation on the turbulent structures, and the interaction between unsteady bubble formation and vortex dynamics. Some studies have been performed using RANS [START_REF] Zhu | Extension of the Schnerr-Sauer model for cryogenic cavitation[END_REF][START_REF] Tseng | Investigations of empirical coefficients of cavitation and turbulence model through steady and unsteady turbulent cavitating flows[END_REF][START_REF] Goncalves | Numerical simulation of cavitating flows with homogeneous models[END_REF] and Large Eddy Simulation (LES) turbulence model coupled with different models for cavitation [START_REF] Roohi | Numerical simulation of cavitation around a two-dimensional hydrofoil using {VOF} method and {LES} turbulence model[END_REF][START_REF] Salvador | Computational study of the cavitation phenomenon and its interaction with the turbulence developed in diesel injector nozzles by large eddy simulation (les)[END_REF][START_REF] Dai | Large-eddy simulations of cavitation in a square surface cavity[END_REF][START_REF] Huang | Large eddy simulation of turbulent vortex-cavitation interactions in transient sheet/cloud cavitating flows[END_REF].

By the way, cavitation model is typically dependent on two types of parameters: first, on some physical parameters, such as for example the number of bubbles, that is not usually well measured; secondly, on some empirical parameters, useful for fitting and calibration procedures with respect to the experimental data. Therefore, model parameters represent an important source of uncertainty. When considering RANS simulation, turbulence model coefficients introduce an additional variability. In fact, the assumption of some empirical coefficients remains someway arbitrary, and sometimes tuned for reducing the distance between experiments and numerical solution. For these reasons, it is crucial to consider this epistemic (since it is due to a lack of knowledge) uncertainty in the problem, thus providing a measure of the variability of the numerical solution, i.e. assessing the quality of the numerical prediction. Then, it is of primary importance in cavitating flows to determine not only a converged numerical solution but also a description of the variability of the solution with respect to the known uncertainties, i.e. providing the statistical moments of the quantities of interest. Moreover, it is not an easy task to well define boundary and initial conditions, because of difficulties encountered in order to control accurately experiments in cavitating flows, thus adding an additional source of uncertainty.

Despite the use of stochastic methods applied to the numerical simulation in fluid mechanics is being more and more diffused, only few studies exist concerning the application of uncertainty quantification tools to cavitating flows, and very few works on the interaction between cavitation and turbulence. In 2000 and 2006 Li et al. [START_REF] Li | Stochastic study of cavitation bubbles near boundary wall[END_REF] proposed a Markov stochastic model to reproduce the random behaviour of cavitation bubble(s) near compliant walls. In 2003, Fariborza et al. [START_REF] Fariborza | Intergranular creep cavitation with time-discrete stochastic nucleation[END_REF] proposed an empirical model for the timediscrete stochastic nucleation of intergranular creep cavities. They assumed nucleation to occur randomly in time, with the temporal behaviour being governed by an inhomogeneous Poisson process. In 2007, Giannadakis et al. [START_REF] Giannadakis | Evaluation of the predictive capability of diesel nozzle cavitation models[END_REF] 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. In 2010, Mishra et al. [START_REF] Mishra | Modeling the coupling of reaction kinetics and hydrodynamics in a collapsing cavity[END_REF] introduced a model of cavitation coupled to deterministic and stochastic chemical reactions of solute chemical species. Wilczynski [START_REF] Wilczynski | Stochastic modeling of cavitation phenomena in turbulent flow[END_REF] and Goel et al. [START_REF] Goel | Surrogate model-based strategy for cryogenic cavitation model validation and sensitivity evaluation[END_REF] performed an uncertainties-based study on some hydrodynamic cavitation model parameters. In particular, Wilczynski [START_REF] Wilczynski | Stochastic modeling of cavitation phenomena in turbulent flow[END_REF] applied a stochastic model to capture the interaction of turbulent pressure field on cavitation nuclei population. Moreover, Goel et al. [START_REF] Goel | Surrogate model-based strategy for cryogenic cavitation model validation and sensitivity evaluation[END_REF] 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 [START_REF] Tseng | Investigations of empirical coefficients of cavitation and turbulence model through steady and unsteady turbulent cavitating flows[END_REF] presents a sensitivity study of the cavitation and turbulence models, where different combinations of empirical coefficients are considered. Recently, Rodio & Congedo [START_REF] Rodio | Robust analysis of cavitating flows in the venturi tube[END_REF] 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 steady cavitating flows by coupling a nonintrusive Polynomial Chaos stochastic method with a cavitating CFD solver.

In particular, a steady solution was obtained, even if in [START_REF] Barre | Experiments and modelling of cavitating ows in venturi: attached sheet cavitation[END_REF] experimental measurements showed an unsteady behavior of the cavitating flow. This was due to the choice of the turbulence model. Note that the present paper includes three novelties with respect to [START_REF] Rodio | Robust analysis of cavitating flows in the venturi tube[END_REF]: i) only unsteady RANS simulations are considered here, thus permitting to observe how the statistics of of the cavitation pocket evolve in time, ii) uncertainties on some parameters of the turbulence models are taken into account, where the interest is to determine whether these uncertainties could be of major/minor importance with respect to the inlet measurements uncertainties, iii) experimental data are directly included in the formulation of an inverse problem in order to calibrate turbulence model parameters.

Then, the first contribution of this paper is to consider the propagation of both some uncertainties on turbulence and cavitation model parameters and inlet conditions by using a non-intrusive stochastic method. In particular, uncertainty are propagated through an unsteady numerical solver simulating a Venturi experimental configuration, thus permitting a detailed comparison between experimental and numerical error bars and a study of the solution variability at different times. Several quantities of interest are considered, in terms of void-ratio, pressure and velocity at different stations (where experimental data are available). Second contribution is related to a Bayesian-based computation, which is applied to the calibration of the turbulent model parameters using directly numerical simulation and the available experimental data on void-ratio and velocity at the different stations. This permits to compute the calibrated Probability Density Function (PDF) for two-key parameters of the turbulence model and their associated statistics.

This paper is organized as follows. Section 2 presents the governing equations for reproducing cavitation and the associated numerical methods for solving the deterministic and stochastic part. In Section 3, the experimental configuration under study is described, including the choices in terms of mesh and characterization of the sources of uncertainty. Section 4 is devoted to present some results. First, the forward propagation problem is introduced considering the whole set of uncertainties. Then, the numerical strategy for calibrating some turbulence model parameters is depicted, together with the associated results. Finally, some conclusions and perspectives are drawn in section 6.

Numerical methods

In this paper, non-intrusive methods for Uncertainty Quantification (UQ) are applied, i.e. the CFD code is considered as a black box, and the statistics are computed on some quantities of interest associated to the output of the simulation.

Note that only 2D computations are performed here. Nevertheless, it was clearly observed in various experiments that the topology of cavitation pocket in Venturis or hydrofoils is fully three-dimensional, involving a transversal component for the re-entrant jet (side-entrant jet). In a previous study [START_REF] Decaix | Investigation of three-dimensional effects on a cavitating venturi flow[END_REF], we have performed 3D simulations of the Venturi geometry using the same cavitation software. The side-wall effects and vortical structures along the corner were clearly observed. 3D effects and spanwise velocity were simulated instantaneously in the Venturi. Yet, about time-averaged profiles in the mid-span plane, 3D and 2D simulations provided quite similar results, also in terms of cavitation pocket and time-averaged velocity profiles. Relying on these results, only 2D simulations are considered in this work, also due to the important increase in the computational cost due to the uncertainty propagation.

CFD for the simulation of cavitation

The code is based on the solving of the one-fluid compressible RANS system coupled with the one-equation turbulence model of Spalart-Allmaras. For low Mach number applications, an inviscid preconditioner is introduced. These equations can be expressed as:

P -1 c ∂w ∂t + div (F c -F v ) = S (1)
w =       ρ ρ V ρE α ρν       ; F c =        ρ V ρ V ⊗ V + P I (ρE + P ) V α V ρν V        ; F v =        0 τ v + τ t (τ v + τ t ). V -Q v -Q t 0 (µ + ρν/σ ν ) grad ν       
where w denotes the conservative variables and the void ratio, F c and F v the convective and viscous flux densities and S the source terms, which concern the void ratio equation and the transport equations. The expression of the preconditioning matrix P c is given in [START_REF] Goncalves | Modelling for isothermal cavitation with a four-equation model[END_REF]. ν is the transported turbulent variable. The total stress tensor τ is evaluated using the Stokes hypothesis, Newton's law and the Boussinesq assumption. The total heat flux vector Q is obtained from the Fourier law involving a turbulent thermal conductivity λ t with the constant Prandtl number hypothesis.

τ = τ v + τ t = (µ + µ t ) ( grad V + ( grad V ) t ) - 2 3 ( div V )I Q = Q v + Q t = -(λ + λ t ) grad T with λ t = µ t C p P rt (2)
In pure phases, the viscosity is assumed to be constant. The mixture viscosity is defined as the arithmetic mean of the liquid and vapour viscosities:

µ(α) = αµ v + (1 -α)µ l (3)
The mixture thermal conductivity λ is also defined as the arithmetic mean of the liquid and vapour values:

λ(α) = α µ v C pv P rv + (1 -α) µ l C p l P r l (4)
The turbulent Prandtl number P rt is set to 1.

To compute the pressure and the temperature, an equation of state (EOS) is necessary to link these thermodynamic quantities to the internal energy and the density. For the pure phases, we used the convex stiffened gas EOS:

P (ρ, e) = (γ -1)ρ(e -q) -γP ∞ (5) P (ρ, T ) = ρ(γ -1)C v T -P ∞ (6)
T (ρ, h) = h -q C p (7)
where γ = C p /C v is the heat capacity ratio, C p and C v are thermal capacities, q the energy of the fluid at a given reference state and P ∞ is a constant reference pressure.

The cavitation model

The cavitation model applied in the present work is based on a void ratio transport equation model:

∂α ∂t + div (α V ) = (K + α) div V + ṁ ρ I (8)
K = ρ l c 2 l -ρ v c 2 v ρ l c 2 l 1-α + ρvc 2 v α ; 1 ρ I = c 2 v α + c 2 l 1-α ρ l c 2 l 1-α + ρvc 2 v α (9)
where ṁ is the mass transfer between phases and ρ I the interfacial density .

By assuming that the mass transfer is proportional to the divergence of the velocity, it is possible to build a family of models in which the mass transfer ṁ is expressed as [START_REF] Goncalves | Numerical study of expansion tube problems: Toward the simulation of cavitation[END_REF]

ṁ = ρ l ρ v ρ l -ρ v 1 - c 2 c 2 wallis div V ( 10
)
where c wallis is the propagation velocity of acoustic waves without mass transfer [START_REF] Wallis | One-dimensional two-phase flow[END_REF]. This speed of sound is expressed as a weighted harmonic mean of speeds of sound of each phase:

1 ρc 2 wallis = α ρ v c 2 v + 1 -α ρ l c 2 l (11)
A first model was built using the speed of sound associated with a sinusoidal barotropic EOS [START_REF] Goncalves | Numerical simulation of cavitating flows with homogeneous models[END_REF][START_REF] Delannoy | Two phase flow approach in unsteady cavitation modelling[END_REF]. Downstream the sheet cavities, it has been observed that the convected void ratio was not enough destroyed. An empirical destruction term was added [START_REF] Goncalves | Modelling for isothermal cavitation with a four-equation model[END_REF]:

ṁ = ρ l ρ v ρ l -ρ v 1 - c 2 c 2 wallis div V -C des ρ v ρ l α Max(0, P -P vap ) 0.5ρ ref U 2 ref ( 12
)
where C des is a tunable parameter.

The turbulence model

The one-equation Spalart-Allmaras model [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] is used in the present study. The transport-equation for ν is:

∂ρν ∂t + div ρν V - 1 σ ν (µ + ρν) grad ν = c b 1 (1 -f t2 ) Sρν + c b 2 σ ν grad ρν. grad ν -c w 1 f w - c b 1 κ 2 f t2 ρ ν2 d 2 (13)
with:

ν t = νf v 1 ; f v 1 = χ 3 χ 3 + c 3 v 1 ; χ = ρν µ S = |rot V | + ν κ 2 d 2 f v 2 ; f v 2 = 1 - χ 1 + χf v 1 ; f t2 = c t3 exp(-c t4 χ 2 ) f w = g 1 + c 6 w 3 g 6 + c 6 w 3 1 6
;

g = r + c w 2 (r 6 -r) ; r = ν Sκ 2 d 2
Constants are:

c b 1 = 0.1355 ; c b 2 = 0.622 ; σ ν = 2 3 ; κ = 0.41 c w 1 = c b 1 κ 2 + (1 + c b 2 )/σ ν ; c w 2 = 0.3 ; c w 3 = 2 c v 1 = 7.1 ; c t3 = 1.1 ; c t4 = 2
Such transport-equation turbulence models always leads to the generation of stable cavities, because strong turbulent eddy viscosity µ t inside the cavity avoids the re-entrant jet formation which plays the major role on the instability of partial sheet cavity. As a remedy to reduce the turbulent viscosity, we use an ad hoc limiter in the mixture area proposed by Reboud [START_REF] Stutz | Two-phase flow structure of sheet cavitation[END_REF][START_REF] Goncalves | Numerical study of unsteady turbulent cavitating flows[END_REF][START_REF] Decaix | Time-dependent simulation of cavitating flow with kℓ turbulence models[END_REF]. This limiter introduces a function f (ρ) in the computation of the turbulent viscosity:

µ t = f (ρ)ν t with f (ρ) = ρ v + (1 -α) N reb (ρ l -ρ v ) ( 14
)
where N reb is a tunable parameter usually fixed to 10.

Wall functions

For the modelling of flow close to the wall, a two-layer wall law approach is used:

u + = y + if y + < 11.13 u + = 1 κ ln y + + 5.25 if y + > 11.13 u + = u U τ ; y + = yU τ ν w ; U 2 τ = τ w ρ w (15)
where κ = 0.41 is the von Karman constant and the subscript 'w' is used for a wall value. We assume that wall functions are similar in a two-phase flow and in a single-phase flow. For unsteady flows, the existence of a wall law is assumed to be valid at each instant. These assumptions have been studied in [START_REF] Goncalves | Wall model and mesh influence study for partial cavities[END_REF] and comparisons were proposed with a thin boundary layer approach.

Numerics

The numerical solving is carried out using an implicit integration based on a finite-volume discretization. For the mean flow, the convective flux density vector on a cell face is computed with the Jameson-Schmidt-Turkel scheme [START_REF] Jameson | Numerical solution of the Euler equations by finite volume methods using Runge-Kutta time stepping schemes[END_REF]. The artificial viscosity includes a second-order dissipation term and a fourth-order dissipation term. The viscous terms are discretized by a second-order space-centered scheme. For the turbulence transport equations, the upwind Roe scheme [START_REF] Roe | Approximate Riemann solvers, parameters vectors, and difference schemes[END_REF] is used for the convective flux. The second-order accuracy is obtained by introducing a flux-limited dissipation. The diffusive flux is discretized by a second-order centered scheme. Time integration is achieved using the dual time stepping approach and a low-cost implicit method consisting in solving, at each time step, a system of equations arising from the linearization of a fully implicit scheme. The derivative with respect to the physical time is discretized by a second-order formula. The turbulence source term needs special treatment [START_REF] Merci | Computational treatment of source terms in two-equation turbulence models[END_REF]. Only the negative part of the source term Jacobian matrix is considered and replaced by its spectral radius. The numerical treatment of boundary conditions is based on the use of the preconditioned characteristic relationships. More details are given in [START_REF] Goncalves | Numerical simulation of cavitating flows with homogeneous models[END_REF][START_REF] Goncalves | Modelling for isothermal cavitation with a four-equation model[END_REF].

Forward Uncertainty Quantification method

The non-intrusive polynomial chaos, as implemented in the NISP (Non Intrusive Spectral Projection) library (see Ref. [START_REF] Congedo | Shape optimization of an airfoil in a bzt flow with multiple-source uncertainties[END_REF]) is used for propagating uncertainties and in order to build the metamodel in the Bayesian approach.

Using this non-intrusive uncertainty quantification tool means that a single deterministic computation is replaced with a whole set of such 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

u(x, t, ξ) = α u α (x, t)Ψ α (ξ), (16)
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 Φ

α i i (ξ i ), that is, Ψ α (ξ) = n ξ i=1 Φ α i i (ξ i ), where each Φ α i i is a polynomial of degree α i , so that the degree of Ψ α is |α| 1 = n ξ i=1 α i .
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

u α (x, t) = u(x, t), Ψ α Ψ α -2 , ( 17
)
where the scalar product is defined by the expectation operator. For practical use, the PC expansions are truncated to degree No

u(x, t, ξ) = |α| 1 ≤No u α (x, t)Ψ α (ξ). ( 18
)
The number of multivariate polynomials Ψ α , that is, the dimension of the expansion basis, is related to the stochastic dimension n ξ and the degree No of polynomials ; it is given by the formula (

n ξ + No)!/(n ξ !No!).
Several approaches can be used to estimate PC coefficients. The approach used in this study is based on quadrature formulae (see Ref. [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF] for details). The PC coefficients are evaluated from a set of abscissas and weights (ξ i , ω i ) by formulae of the form

u α (x, t) = Ψ α -2 n i=1 u(x, t, ξ i )Ψ α (ξ i )ω i , (19)
where n refers to the number of quadrature points, depending on the quadrature formula. 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 [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF]. In particular, ANOVA relies on the decomposition of the variance according to the contribution of each uncertainty. This permits to identify the contribution of a given stochastic parameter to the total variance of an output quantity, thus computing a hierarchy of the most predominant parameters in terms of the output variability. This can be done obviously for each quantity of interest. In a PC framework, It can be easily computed by using some interesting properties of the previous development. For more details, Ref. [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF] is strongly recommended.

Experimental configuration

The Venturi was tested in the cavitation tunnel of the CREMHyG (Centre d'Essais de Machines Hydrauliques de Grenoble). It is characterized by a divergence angle of 4 • , illustrated in Figure 1. The edge forming the throat of the Venturi is used to fix the separation point of the cavitation cavity. This geometry is equipped with five probing holes to allow various measurements such as the local void ratio, instantaneous local speed and wall pressure (Figure 1).

The selected operating point is characterized by the following physical parameters [START_REF] Barre | Experiments and modelling of cavitating ows in venturi: attached sheet cavitation[END_REF]: U inlet = 10.8 m/s, the inlet velocity

σ inlet = P inlet -P vap 0.5ρU 2 inlet ≃ 0.
55, the cavitation parameter in the inlet section

T ref ≃ 293K, the reference temperature L ref =252 mm, the reference length Re L ref = U inlet L ref ν = 2.7 10 6
, the Reynolds number With these parameters, a cavity length L ranging from 70 mm to 85 mm was obtained. The experimental views for this geometry show a quasi stable cavity behaviour. The attached cavity length corresponding to the end of the re-entrant jet is around 30-35 mm. For this geometry, no periodic cycles with large shedding were observed.

Mesh and numerical parameters

The question of grid convergence in turbulent RANS simulations involving a separation of the boundary layer is difficult to lead because of the problem of wall modelling and turbulence transport-equation stiffness in fine grids. For small first-cells y+ values, the equations for k an ε becomes very stiff and lead to divergence. Large questions are also linked to the wall modelling, wall functions and near-wall turbulence damping functions. These questions become even more difficult for two-phase turbulent boundary layer with separation and cavitation. The existence of a logarithmic law or an universal velocity profile is not attested in such flows. Moreover, two-phase assumptions (such as negligible tension surface effects) are no more valid in fine grids where the size cells can be smaller than 1 micron. In consequence the grid convergence in such cavitating flows raised lots of problems. In a previous work [START_REF] Goncalves | Wall model and mesh influence study for partial cavities[END_REF], we have investigated the mesh influence for the Venturi geometry in the near-wall area when a cavitation pocket develops and interacts with the turbulent boundary layer. We have considered five different meshes, with a final choice consisting in a H-type grid, which contains 251 nodes in the flow direction and 62 nodes in the orthogonal direction. A special contraction of the mesh is applied in the main flow direction just after the throat to better simulate the two-phase flow area. The y + values of the mesh, at the center of the first cell, vary between 12 and 27 for a non cavitating computation.

Unsteady computations are performed with the dual time stepping method and are started from the non cavitating numerical solution. The numerical parameters are:

• the dimensionless time step, ∆t * = ∆tU inlet L ref = 4.88 10 -3
• sub-iterations of the dual time stepping method, 50

• the CFL number, 0.2.

Sources of uncertainty

Two sources of uncertainties are considered in this paper. The first one is related to the whole set of data, coming from several experimental campaigns (see for example [START_REF] Barre | Experiments and modelling of cavitating ows in venturi: attached sheet cavitation[END_REF]) : void ratio and velocity at the five stations described in the previous section, pressure fluctuations at the wall averaged in time, inlet conditions measurements in terms of mass flow rate. For each data, an estimation of the uncertainty has been done in [START_REF] Barre | Experiments and modelling of cavitating ows in venturi: attached sheet cavitation[END_REF]. The velocity is evaluated as the most probable value and the void ratio is obtained from the signal of the double optical probe using a post-processing algorithm. The relative uncertainty on the void ratio measurement were estimated at roughly 15% [START_REF] Barre | Experiments and modelling of cavitating ows in venturi: attached sheet cavitation[END_REF]. Concerning the velocity, uncertainty is assumed to vary according to the law provided in [START_REF] Barre | Experiments and modelling of cavitating ows in venturi: attached sheet cavitation[END_REF]. Pressure at the wall averaged in time is estimated at 0.027% of variation. Concerning the inlet pressure, a maximal variation of 0.05 Bar is considered. Note that inlet/outlet pressure conditions are imposed in the numerical code. The measurement on the inlet pressure P inlet is taking into account as the uncertain inlet condition, and a fixed outlet pressure is imposed, set to 70680 Pa.

The second source of uncertainty consists in some epistemic parameters, used in the cavitation and turbulence model. In particular, we consider as uncertain the Reboud eddy-viscosity limiter N reb , and the destruction parameter C des . For these epistemic uncertainties, relying on the existing literature, the following ranges are considered: intervals [START_REF] Brennen | Cavitation and Bubble Dynamic[END_REF][START_REF] Mishra | Modeling the coupling of reaction kinetics and hydrodynamics in a collapsing cavity[END_REF] and [0, 10] for N reb and C des , respectively. For turbulence model uncertainties, the present analysis is focused on the two parameters which were added ad-hoc in [START_REF] Stutz | Two-phase flow structure of sheet cavitation[END_REF] for considering the cavitation effect (see Eqs. 12 and 14). Note that a complete study should include all the other empirical coefficients of the Spalart-Allmaras turbulence model, even though this could strongly increase the computational cost. This increased cost is related to both compute a complete ANOVA analysis and the difficulty to build a converged response surface if a large number of inputs are considered, which could impact the feasibility of the Bayesian calibration. Nevertheless, the idea is to calibrate the turbulent parameters which are specific to cavitation and are not well-studied in literature. For this reason and because of the potential strong increase of the computational cost, only uncertainties on N reb and C des are considered.

Concerning the PDF, we have used systematically uniform pdf, seeing that an accurate estimation of probability density function for the physical measures, used as input parameters in the numerical simulation, is not available. This choice represents a robust strategy in order to analyze uncertainty propagation of physical uncertainties, since that all the solutions have the same probability. As a consequence, the worst condition is assumed, thus ensuring a robust estimation of the solution variability. Concerning model uncertainty, this epistemic variable is treated again as a uniform pdf, that is one of the possible options when considering this kind of uncertainty.

Two types of analysis are performed in this work. In the first one, epistemic and inlet condition uncertainties are propagated through the numerical CFD solver (the ensemble of uncertainties, i.e. on P inlet , N reb and C des , are resumed in Table 1). The output statistics varying in time, with their numerical error bars, are compared with the experimental measurements with their associated tolerances. In the second analysis, some experimental data, i.e. on the void ratio and the velocities at different stations, are directly used in a Bayesian framework for calibrating the two epistemic uncertainty, on N reb and C des . This permits to compute probability density functions and to statistically characterize these two parameters.

Results

The cavitation numerical code has been already verified, see for more details [START_REF] Goncalves | Modelling for isothermal cavitation with a four-equation model[END_REF]. In the following, a reference deterministic solutions is used for comparison, which has been obtained by assuming N reb = 8, C des = 0.

Forward propagation of the three sources of uncertainties

Uncertainties on N reb , C des and inlet conditions are considered together. Convergence in the stochastic space for mean and variance has been verified by increasing the polynomial order. This study is not reported here for brevity. Each solution described in this section has been obtained with a polynomial order of 6. In the following section, a detailed analysis on the convergence of the response surface coming from the polynomial chaos is presented, since it is directly used during the Bayesian approach. In this section, mean and variance for different quantities of interest are computed for each time-step, thus permitting to estimate the variation in time of the cavitating flow. This is a fundamental difference with respect to [START_REF] Rodio | Robust analysis of cavitating flows in the venturi tube[END_REF], where only steady computations have been considered. Then, in this case the simulation of the flow description was not sufficiently accurate to describe the experimental flow behavior observed in [START_REF] Barre | Experiments and modelling of cavitating ows in venturi: attached sheet cavitation[END_REF]. In this section, along with the description of the current results, a systematic comparison with the steady results of [START_REF] Rodio | Robust analysis of cavitating flows in the venturi tube[END_REF] is done in order to highlight the interest of an unsteady study under uncertainties.

First, mean (in the stochastic sense, i.e. with respect to the considered uncertainties) for the pressure is computed and compared to the reference solution (see Figure 2) at the final time-step. A low-pressure area is clearly observed downstream of the throat (abscissa x = 0 m) corresponding to the sheet cavity. The vapor pressure value is set to its physical value at the reference temperature: P vap = 2340 Pa. For the deterministic solution, the low-pressure area is more extended in comparison with the stochastic mean solution. At the closure part of the cavity (around abscissa 0.08-0.1 m), small cavitation clouds are shed (see the next figure) that explained the low values of the pressure field. This behaviour is not observed for the stochastic mean solution. Downstream, in the recompression area, the condensation process is well calculated. Secondly, the void ratio is then considered and the solution computed for four different times (times equal to t=1.14, 1.5, 1.9 and 2.28 s). Also in this case, solution is compared with the deterministic reference one (see Figure 3 on the right). The sheet cavity is composed by an attached part fixed to the throat with high values of void ratio (up to 0.98) following by a fluctuating recirculation area involving small cavitation clouds shedding. The re-entrant jet phenomenon is not enough developed to clearly break off the pocket and to generate large clouds shedding. This dynamics behaviour simulated by the reference solution is in close agreement with the experimental visualization. As regard to the stochastic mean solution, the sheet behaviour is quite different, except at the initial time t = 1.14 s. The evolution of the mean void ratio put in evidence the extension of the mean cavity length, which reaches the abscissa x = 0.2 m at the final time. Moreover, with such a cavity length, we expect that the low pressure mean values follows a similar shape, that is not the case (see the previous figure). These two quantities seem to be not correlated contrary to the deterministic quantities. The different stochastic behaviour with respect to the deterministic one is motivated by the strong influence that uncertainties have on the development of the sheet for increasing times.

Let us now focus on the variance of the void ratio σ 2 (α), which is plotted in Figure 4 at the same four times. As observed previously on the mean void ratio, the variance is very weak at the initial time and increases with the time marching. This means that higher values of void-ratio are obtained for different combinations of the considered uncertainties (looking mainly at the means), but also that more different values of void-ratio are computed. This makes the prediction on the void-ratio questionable since higher means are associated to higher variance, thus yielding a large variability of the predicted solutions.

Decomposing the variance by means of an ANOVA analysis, permits to compute the contribution of each source of uncertainty to the global variance. In this case, two times are considered, for evaluating the evolution with respect to the sensitivity indices: t = 1.5 s and the final time t = 2.28 s. Results are presented in Figures 5 and6 for the two times, respectively. The same qualitative results are found in both cases. For both model parameters Cdes and N Reb , values of sensitivity indices are very low in the two-phase mixture area and the maximum value is around 0.4%. Then, the influence of the epistemic uncertainties associated to the model are someway negligible in terms of variance production. On the other hand, the pressure inlet sensitivity index is clearly dominant with values 50 times higher. This means that experimental uncertainties on inlet boundary conditions are predominant with respect to the model-uncertainty.

Note that the comparison with results known in literature [START_REF] Rodio | Robust analysis of cavitating flows in the venturi tube[END_REF] since these last ones deals with steady computations. Nevertheless, a strong similarity can be observed between the mean and the variance of the void ratio reported in [START_REF] Rodio | Robust analysis of cavitating flows in the venturi tube[END_REF], with the solution computed here at a time between 1.14 and 1.5 s. Concerning the ANOVA contributions, similar conclusions from a qualitative point of view were done in [START_REF] Rodio | Robust analysis of cavitating flows in the venturi tube[END_REF] based on steady simulations for the same geometry. In fact, even if a quantitative comparison is not possible because of the different turbulence and cavitation models used (note that in [START_REF] Rodio | Robust analysis of cavitating flows in the venturi tube[END_REF] no uncertainty on turbulence model parameters are considered), inlet uncertainties are predominant with respect to the model ones. As a consequence, numerical prediction could be substantially improved by reducing the aleatory uncertainty associated to the measurements at the inlet conditions. Another quantity of interest is the RMS wall pressure fluctuations, i.e. P ′ rms , reported in Figure 7. The pressure fluctuation is divided by the timeaveraged pressure, P av . We can observe that at the end of the sheet cavity, close to x=0.19 m, the experimental pressure fluctuations increase to a value of 0.2, but, at the downstream of the cavity, they decrease to zero. The mean of the numerical pressure fluctuations shows a perfect agreement with the experimental data in the region of the sheet cavity. However, a numerical over-estimation is evident in the downstream region. This same trend in terms of fluctuations, when a Spalart-Allmaras model is used, has been observed in previous studies (see for example [START_REF] Charriere | A comparative study of cavitation models in a venturi flow[END_REF]), where the influence of several turbulence models has been assessed. This is a real issue when SA turbulence model is considered. Note also that numerical error bars due to the uncertainties are very small. Because of this lack of prediction due to the structural form of the model, experimental data on pressure could not be considered in the calibration of model parameters. This assumption is also explained in Section 4.2. The sensitivity indices for the pressure fluctuations are also plotted in Figure 7 in order to determine the contribution of each source of uncertainty to the global variance for this quantity. It clearly shows that the N reb turbulence parameter contribution is predominant with respect to the others, thus confirming that the importance of turbulence model in the description of the RMS wall pressure fluctuations. In Figure 8, the mean cavitation number, σ = (P inlet -P vap )/0.5ρu 2 inlet , is plotted versus the physical time. Usually, the cavitation number is a good indicator about the onset and the evolution of the cavitation regime. In particular, a decrease of σ corresponds to an increase of the cavitation sheet. In this case, a non-negligible increase of the mean of σ can be observed, even if the mean of the cavitation sheet grows in time (see Figure 4), starting from a time of 1.48 s. Considering that the cavitation number is dependent mostly from the inlet pressure, P inlet , the behaviour of the σ seems to indicate that there is not a correlation between the means of size of the cavitation sheet and σ. Note anyway the enormous variations in terms of coefficient of variation, i.e. the ratio of the standard deviation to the mean, that can attain also the 20%, that indicate this quantity is not predicted in a robust way. Let us now compare the numerical error bars of the void ratio and of the velocity at each station (described in Section 3) with the available experimental measurements (see Figures 9 and10, respectively), at the final time-step. Concerning the void ratio, note first that, quantitatively, means (numerical and experimental ones) are quite distant, then it is hard to get interesting information looking only at the means. Moreover, a good estimation could not be done by performing a simple deterministic simulation without considering the strong variability due to the uncertainties. Note also that looking at the numerical and experimental error bars, they are crossing each other. This means that the numerical solution is someway a "conservative" approximation of the considered quantity of interest. In fact, experimental error bars in several points lie inside the numerical ones. Finally, note that numerical error bars are very large. This means that for some points, nearly every possible value of α is possible. It seems then of prominent importance to improve the numerical prediction, by trying to reduce the sources of uncertainty affecting the numerical simulation.

t (s)

Though only steady computations are provided in [START_REF] Rodio | Robust analysis of cavitating flows in the venturi tube[END_REF], void ratio trends in stations 1 and 5 are very similar in terms of mean and variance. On the contrary, errors at station 3 are different, and larger error bars are observed here, probably due to the use of a different turbulence model.

Let us consider now the velocity, represented in Figure 10 at the different stations. At the first station, an inconsistency is observed between the numerical and experimental solution. Looking at the experimental values, it seems questionable that velocity at y = 0.0012 m could decelerate, but additional experimental campaign should be performed for verifying this point and thus validating the numerical solution. Nevertheless, numerical error bars at station 1 are small, even smaller than the experimental bars. From station 2 to 5, numerical prediction is quite accurate, except getting close to the wall. Note anyway that in these cases, numerical error bars become larger. Globally, velocities are not captured with the same notion of conservative bars, mentioned for the void-ratio, which is maybe due to a deficiency in the model. Note that a recirculation region exists (associated to the values of negative velocity) at the fourth and fifth stations, which is only partially well-captured. Comparing with the steady results obtained in [START_REF] Rodio | Robust analysis of cavitating flows in the venturi tube[END_REF], note that velocity profiles are very similar, except in station 5, where the turbulence model used in [START_REF] Rodio | Robust analysis of cavitating flows in the venturi tube[END_REF] seem to provide better results in terms of comparison with experimental data.

Parameter Calibration with respect to the experimental data.

The output of Bayesian inference is not a single value for the model parameters, but a posterior probability distribution that summarizes all available information about parameters. From this distribution, one can calculate means, modes, and high-order moments, compute marginal distributions, or make additional predictions by averaging over the posterior.

Let F denote the forward mathematics model defined as follows : d = F (m), which yields predictions of the void-ratio and velocity at different stations d = (α, u) as a function of the Reboud eddy-viscosity limiter N reb , and the destruction parameter C des , i.e. m = (N reb , C des ). In the Bayesian setting, m are random variables and we use Bayes' rule to define a posterior probability density for the model parameters m, given n observations of the data {d 1 , . . . , d n } :

p(m|d 1 , . . . , d n ) = p(d 1 , . . . , d n |m)p m (m) p(d 1 , . . . , d n |m)p m (m)dm . ( 20
)
Prior probability p m (m) represents the degree of belief about possible values of m = (N reb , C des ) before observing any data ; N reb and C des are a priori assumed to follow uniform distributions, with minima and maxima defined in Section 3.2.

Data then enters the formulation through the likelihood or joint density of the observations given m, namely p(d 1 , . . . , d m |m). A common model assumes independent observations so that independent additive errors account for the deviation between predicted and observed values of d :

d j = F (m) + η j , j = 1, . . . , n. (21)
Because α and u can be considered as independent, a typical assumption is that errors are realizations of a Gaussian random variable η j ∼ N (0, Γ), where Γ = diag(σ 2 α , σ 2 u ), σ α and σ u encompassing measurement errors. In that case, d j |m ∼ N (F (m), Γ), and the likelihood is

p(d 1 , . . . , d m |m) = n j=1 p d j (d j |m) = n j=1 p η (d j -F (m), Γ), (22)
with p η the Gaussian density probability of N (0, Γ). Since in general measurement errors are not known with exactness, one can consider σ st = (σ α , σ u ) as hyperparameters of the Bayesian setting that needs to be inferred, with noninformative uniform a priori on σ α and σ u . In this case, measurement errors are known, and are defined in Section 3.2. Markov Chain Monte Carlo (MCMC) encompasses a broad class of methods that simulate drawing samples from the normalized posterior [START_REF] Gilks | Markov Chain Monte Carlo in Practice[END_REF]:

p(m|d 1 , . . . , d n ) ∝ p(d 1 , . . . , d n |m)p m (m), (23)
thus avoiding complex numerical integrations in high dimensions to form the posterior distribution. In this work, the Metropolis-Hastings algorithm with single-site updating and Gaussian proposal density to draw samples of p(m|d 1 , . . . , d n ) and process is used (see [START_REF] Tryoen | Bayesian-based method with metamodels for rebuilding freestream conditions in atmospheric entry flows[END_REF] for more details).

A strong issue is that we can not afford to call the CFD cavitating code two times for each iteration of the Markov Chain, since the code is timeconsuming (about one hour per simulation) and one needs some thousands of iterations to produce a good sample of the posterior. To tackle this issue, one can rely on a metamodel, which gives an approximation of the outputs of the Cavitation Code as a function of its inputs. Note that the number of samples for building an accurate metamodel is orders of magnitude less than the number of iterations of MCMC. Metamodels based on intrusive and non-intrusive stochastic spectral methods have already been proposed in the context of Bayesian inference [START_REF] Marzouk | Stochastic spectral methods for efficient Bayesian solution of inverse problems[END_REF][START_REF] Marzouk | A Stochastic Collocation Approach to Bayesian Inference of Inverse Problems[END_REF], with PC expansions as presented in Section 2.2.

PC expansion metamodel functions of the form [START_REF] Rodio | Robust analysis of cavitating flows in the venturi tube[END_REF] are obtained, providing an approximation of the response (α, u) as a function of (m), in particular for each physical location at the five measurements stations. In order to build the metamodel, a full tensorized quadrature formula of level 6, requiring 343 resolutions of the deterministic code is performed. Then, different polynomial order are considered, from 3 to 6. In Figures 11 and12, mean and variance of the void-ratio for each polynomial order at stations 1 and 3 are reported, respectively. As it can be observed, a very good convergence is attained at order 6. Similar trend are observed also for the different stations not reported here for brevity. Moreover, a sixth polynomial order permits a very good convergence also for the velocity. For choosing the polynomial order, two assessments are performed. First, approximated response surfaces of α and u at different stations are represented as a function of N reb , C des and P inlet and compared with the outputs provided by the resolutions of the cavitation code. Error norms on the void-ratio and velocity are below the 0.07%. Secondly, probability density functions for the different orders are compared in order to estimate if artificial oscillations are introduced when using a higher polynomial order. Finally, surrogate models retained for the inverse problem are built using a polynomial order No = 6 for both the void-ratio and the velocity.

As a consequence, the stochastic inverse problem is solved by considering the measurements on the void-ratio and the velocity at the five stations (described in Section 3). Note that measurements on the pressure are not used in the inverse problem, because of the behavior of the cavitation model for simulating pressure signals (see Section 4.1). For details concerning the MCMC algorithm, refer to [START_REF] Tryoen | Bayesian-based method with metamodels for rebuilding freestream conditions in atmospheric entry flows[END_REF].

Applying the MCMC algorithm allows solving the problem described in Eq. 20, focused on the calibration of the two empirical parameters of the cavitation and turbulence model, i.e. N reb and C des . Results showing the chain position over 30000 iterations are reported in Figure 13, plotted in two dimensions and separately for N reb and C des . Visual inspection suggests that the chain mixes well and that it moves in the band of where the probability is nonzero.

Once the algorithm is converged, the Probability Density Functions (PDF) can be computed. In particular, Figure 14 represents the computed PDF for N reb and C des . From the MCMC sample, the posterior distribution can be also summarized in terms of means, standard deviations, and marginal distributions [START_REF] Gilks | Markov Chain Monte Carlo in Practice[END_REF]. In order to approximate them, it is important to drop the first (b -1) iterations of the MCMC sample, where b is commonly called the "burn-in" time, necessary for the chain to reach a good behavior. In Table 2, some quantities, i.e. the mean, the coefficient of variation and the 95% confidence intervals, are reported. Note that the PDF for N reb is quasi-gaussian, with a coefficient of variation of nearly the 1.3%, thus providing some good insights for the optimal choice of this parameter. On the other hand, the PDF for C des displays a quasi multi-modal behaviour. Note also that in this case, the coefficient of variation (and the 95% confidence intervals, is very large (nearly 28%), thus confirming the widespread distribution of the PDF. Finally, concerning the parameter C des , it is questionable to indicate an optimal value to use, because of its large variability.

Mean

Coefficient of variation 95%

Conclusions

This paper is focused on the uncertainty assessment of cavitating flows in terms of cavitation and turbulence models and inlet conditions. In particular, uncertainties are propagated through an unsteady numerical solver simulating a Venturi experimental configuration, thus permitting a detailed comparison between experimental and numerical error bars and a study of the solution variability at different times. Furthermore, a Bayesian-based approach has been applied for the calibration of model parameters using fluids numerical simulation and the available experimental data on void-ratio and velocity at the different stations. This permits to compute the calibrated Probability Density Function (PDF) for two-key parameters of models and their associated statistics.

This study has clearly shown how the uncertainty characterization and the computation of the stochastic variability is of fundamental importance for the assessment of the proposed physical model, and for determining which quantities of interest could be harder to predict. Being able to estimate the predictive character of the simulation is even more important in complex flows, where the choice of modelling could be very challenging, also due to scarce experimental data.

This approach has been applied on the simulation of a Venturi configuration. Several quantities of interest have been considered, in terms of voidratio, pressure and velocity at different stations (where experimental data are available), by analyzing also the evolution in time. First, it has been observed that the mean cavity length is someway decorrelated to the evolution of the mean pressure, which is contrary with respect to the deterministic quantities. Secondly, by decomposing the variance by means of an ANOVA analysis, the influence of the epistemic uncertainties associated to the model has shown to be negligible in terms of contribution to the global variance. As a consequence, experimental uncertainties on inlet boundary conditions are predominant with respect to the model-uncertainty. Furthermore, not a good agreement has been observed for the RMS wall pressure fluctuations, except in the region of the sheet cavity. By comparing the numerical error bars of the void ratio and of the velocity at each station, it can be observed that the numerical solution can not be considered as predictive, since very large error bars are computed. Nevertheless, the numerical solution has provided a conservative approximation of the void-ratio, since experimental error bars in several points lie inside the numerical ones. Concerning the velocity, a general agreement is observed, but with some inconsistencies between the numerical and experimental solutions, probably due to the turbulence model. Finally, a bayesian Calibration of turbulent model parameters has been performed using available experimental data on the void-ratio and the velocity for each station. Results has shown that the obtained PDF for N reb is quasi-gaussian, with a small coefficient of variation, thus providing a good indication for the choice of this parameter. On the other side, the PDF for C des displays a quasi multi-modal behaviour. Due to the large confidence interval, it is questionable to indicate an optimal value to use for this parameter.

Future studies will be directed towards the investigation of different turbulence models with their associated uncertainties. The bayesian calibration should be applied for different families of models for providing an estimation of the optimal model to use for this kind of flows. An important aspect that has not been considered in this paper is also the computation of the numerical error associated to the mesh. Here, only a convergence mesh study has been performed on the basic deterministic numerical simulation. Nevertheless, the error coming from the mesh, could be non-negligible, for instance when cavitation sheets increase in time. For this reason, a future aim is to implement a specific strategy, for example the one presented in [START_REF] Ea | A procedure for the estimation of the numerical uncertainty of cfd calculations based on grid refinement studies[END_REF], for considering the mesh effect as an additional source of uncertainty. Another relevant point coming from the analysis presented in this study could concern new experimental investigation for trying to reduce the uncertainties on the inlet conditions.

Acknowledgement

Figure 1 :

 1 Figure 1: Experiment.

Figure 2 :

 2 Figure 2: Reference (top) and mean pressure (bottom) at the final time-step.

Figure 3 :

 3 Figure 3: Mean (left) and reference void ratio (right) at four different times.

Figure 4 :

 4 Figure 4: Variance of the void ratio at four different times.

Figure 6 :

 6 Figure 6: Sensitivity indices of the void ration at a time t = 2.28 s.

Figure 7 :

 7 Figure 7: RMS wall pressure fluctuations along the profile, comparison between experimental data and numerical solutions (on the left), ANOVA analysis (on the right).

Figure 8 :

 8 Figure 8: Time evolution of σ cav .

5 Figure 9 : 5 Figure 10 :Figure 11 :Figure 12 :

 595101112 Figure 9: Comparison of numerical and experimental error bars for α at the different stations described in Section 3 at the final time-step.

Figure 13 :Figure 14 :

 1314 Figure 13: MCMC chain over 30,000 iterations, plotted in two dimensions for N reb and C des , respectively.

Table 1 :

 1 Deterministic values and maximal variation for model and operating conditions uncertainties

		Deterministic values Max variations
	P inlet	36000 (Pa)	±5000 (Pa)
	N reb	8.0	2.0 -15.0
	C des	0.0	0.0 -10.0

Table 2 :

 2 Some statistics computed from the PDF for N reb and C des .

The authors thank the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-11-MONU-0008.