Bayesian Inference of Cavitation Model Coefficients and Uncertainty Quantification of a Venturi Flow Simulation

In the present work, uncertainty quantification of a venturi tube simulation with the cavitating flow is conducted based on Bayesian inference and point-collocation nonintrusive polynomial chaos (PC-NIPC). A Zwart–Gerber–Belamri (ZGB) cavitation model and RNG k-ε turbulence model are adopted to simulate the cavitating flow in the venturi tube using ANSYS Fluent, and the simulation results, with void fractions and velocity profiles, are validated with experimental data. A grid convergence index (GCI) based on the SLS-GCI method is investigated for the cavitation area, and the uncertainty error (UG) is estimated as 1.12 × 10−5. First, for uncertainty quantification of the venturi flow simulation, the ZGB cavitation model coefficients are calibrated with an experimental void fraction as observation data, and posterior distributions of the four model coefficients are obtained using MCMC. Second, based on the calibrated model coefficients, the forward problem with two random inputs, an inlet velocity, and wall roughness, is conducted using PC-NIPC for the surrogate model. The quantities of interest are set to the cavitation area and the profile of the velocity and void fraction. It is confirmed that the wall roughness with a Sobol index of 0.72 has a more significant effect on the uncertainty of the cavitating flow simulation than the inlet velocity of 0.52.


Introduction
Recently, uncertainty quantification (UQ) in computational fluid dynamics (CFD) has become an important issue in the improvement of computational performance and the reduction in computational cost in many engineering problems, such as hypersonic flow simulation for aero vehicles [1], heat and fluid flow simulations in nuclear reactors [2], and blood flow simulation in the vessels [3], etc.
Uncertainty quantification methodology can be divided into two types, depending on the presence or absence of a feedback process of statistical random inputs-the forward problem and the inverse problem [4]. The forward problem sets up the probability distribution of random inputs, calculates the amount of propagation of the inputs to the output, and calculates the quantities of interest (QoI) based on the selected UQ method. However, in the reverse problem, Bayesian inference is applied to estimate the statistics of random inputs, called a posterior distribution, using observation data with experimental or synthetic data.
In the case of the forward problem, Shaefer et al. [5] adopted the point-collocation nonintrusive polynomial chaos (PC-NIPC) method for aerodynamic flow simulation with turbulence models. They [5] investigated the effect of sampling a number of random inputs depending on the polynomial order and the probability distribution of aerodynamic performance according to the statistics of turbulence model coefficients. The model coefficients of the Spalart-Allmaras (SA) turbulence model and k-ω family models were assumed to be uniform distributions with reasonable intervals, and the geometry of the RAE 2822 airfoil and an axisymmetric bump were considered. The coefficients that contribute most to uncertainty in the results were presented according to each model and the adopted geometry.
Cheung et al. [4] proposed the methodology and process of Bayesian uncertainty quantification for the calibration and prediction of the model coefficients and combined them with a CFD simulation. They [4] investigated the effect of stochastic model classes with no uncertainty, independent Gaussian uncertainty, and correlated Gaussian uncertainty. The coefficients of the SA turbulence model were calibrated from the simulation of a flat plate boundary layer and its experimental data. They used a deterministic boundary layer solver directly instead of the surrogate model for the adaptive multilevel algorithm based on Markov chains. Their results emphasized the importance of model inadequacy and observation error in the uncertainty model, and the model class with correlated Gaussian uncertainty offered the best probabilistic results.
OCED Nuclear Energy Agency (NEA) published a review report on the UQ methodologies and the applicability of proposed UQ in CFD for nuclear safety and the regulation of equipment and components [6]. They divided the category of the possible application of single-phase CFD for Nuclear Reactor Safety (NRS) into domains and proposed a degree of maturity to each UQ method and proper choice depending on the subdomain. In Chapter 10 of the report [6], one case of uncertainty quantification for mixing flow CFD results was introduced. This forward problem adopted the nonintrusive generalized polynomial chaos expansion (GPCE) method and estimated the uncertainty of the predicted velocity, concentration, and turbulent quantities through comparison with experimental data in the Paul Scherer Institute GEMIX facility.
In the present work, the venturi tube, which is one of the important components in the field of in-service testing during reactor operation, is considered for Bayesian uncertainty quantification combining the CFD method.
The venturi tube includes cavitating flow where the vapor bubble is formed due to the phase change by a local pressure drop. In addition, the surface roughness may affect the venturi tube's inherent function and operability of the safety-related system. However, the accurate prediction of the cavitating flow is still a challenging issue because of the model's inadequacy in CFD and the measurement difficulties in experiments. This implies that the uncertainty of cavitating flow simulation inside a venturi tube installed at the in-service testing-related systems should be investigated with respect to the adopted cavitation model and operating conditions.
Barre et al. [7] conducted both the experiments and the numerical simulations on the cavitating flow with attached sheet cavitation in a venturi geometry. They proposed a new double optical probe measurement methodology and vapor detection algorithm to measure the void fraction and velocity profiles in the venturi tube. For the computational simulation, a barotropic approach and Yang and Shih k-ε turbulence model [8] were adopted. They showed that the simulation results are in good agreement with experimental data in a stable zone of cavitation but that there were some discrepancies in the flow velocity and void fraction profiles in the rear part of the cavitation sheet. Additionally, the sensitivity analysis of AMIN, the parameter which links the fluid density to the local pressure in the barotropic model, is conducted for better prediction of cavitating flow features [7].
Rodio and Congedo [9] applied a nonintrusive stochastic chaos method to the CFD solver to analyze the uncertainty affecting the prediction of cavitation flow in a venturi tube. The Schnerr and Sauer cavitation model [10] and k-ε turbulence model were employed. For uncertainty quantification, one model coefficient (η) and two flow conditions (inlet pressure and velocity) were considered as random variables. They found that the inlet boundary conditions are dominant contributors to the model uncertainty. In addition, a simple algorithm to obtain the optimized parameter of epistemic uncertainty was proposed. Goel et al. [11] performed the model validation in a transport-based cryogenic cavitation model for a simple 2D hydrofoil and studied the cavitation model parameter and uncertainties with respect to temperature using surrogate-based global sensitivity analysis. They found that the model parameter associated with the evaporation was more dominant than the parameter of material uncertainty. The optimized model parameters of the transport-based cavitation model were presented in the cryogenic cavitation flow.
Ge et al. [12][13][14] investigated the thermal effect on hydraulic cavitation dynamics through experiments. They confirmed that the thermal effects showed a significant role and inhibited the growth of the sheet cavity as the temperature increased. Moreover, it was found that the re-entrant jet triggered the detachment of the attached cavity and unsteadiness.
In the current work, a CFD simulation for the venturi tube flow is combined with the Bayesian inference methodology to quantify the uncertainty of simulation results. First, the cavitation model coefficients are calibrated based on Bayesian inference with the experimental data from observations. Next, using the calibrated model coefficients from the previous section, the uncertainties of the simulation results are quantified using the PC-NIPC method for the forward problem. In this case, two operating conditions of the inlet velocity and the wall roughness at the venturi tube boundary are considered uncertain parameters, and the QoI are the cavitation area and the profiles of void fraction and velocity. In this regard, one of the present co-author [15,16] successfully predicted the location of flow leakage caused by the cavitation erosion and found that the averaged vapor volume fraction tended to increase as the surface roughness decreased in the multistage orifice.
This paper is organized as follows. Section 2 presents the governing equation for the simulation of the cavitation flow in the venturi tube and the overall procedure of Bayesian inference. The PC-NIPC model for the construction of the surrogate model is explained briefly. Then, the procedure to calculate the likelihood function, including model inadequacy and observation error, is presented. In Section 3, the venturi geometry and operating conditions are described, and the deterministic solution is validated through comparison with experimental data. To quantify the uncertainty of the grid resolution, the grid convergence index (GCI) is analyzed with respect to coarse, medium, and fine meshes. Section 4 presents the results of the Bayesian inference of ZGB model coefficients. In Section 5, the two operating conditions of inlet velocity and wall roughness are considered, and the uncertainty of the cavitation area is quantified for the forward problem.

Cavitation and Turbulence Model
2D Reynolds-averaged Navier-Stokes equations are considered using the commercial CFD code ANSYS FLUENT v.19.1 [17]. To predict the cavitation phenomenon, a mixture multiphase model and Zwart-Gerber-Belamri (ZGB) cavitation model [18] are adopted. The turbulence model and wall modeling in the present work are selected based on the previous research of Barre et al. [7] and Rodio and Congedo [9]. The renormalization group (RNG) k-ε model is adopted for the turbulence model. The near-wall region is modeled by a standard wall function. The SIMPLE algorithm [19] is used for the pressure-velocity coupling, and a second-order upwind scheme is used for spatial discretization. For a better convergence rate, the pseudo-transient algorithm is adopted.
The governing equations of the ZGB model are as follows: where f v is the vapor mass fraction and Γ is the diffusion coefficient. R e and R c represent the production and destruction terms in cavitation, respectively, with 4 model coefficients as shown in Table 1. R B and α nuc are the bubble radius and nucleation site volume fractions, respectively. F vap and F cond are evaporation and condensation coefficients which are multiplied in the production and destruction terms.

Grid Convergence Index (GCI)
A grid is generated to discretize the calculation domain in CFD. The number and quality of the grid elements will affect the accuracy and the errors in the calculation results. Therefore, an appropriate level of the grid should be considered, and GCI is used as a method to quantify the uncertainty according to the grid qualities. The GCI is an indicator of how the solution will converge with a further grid refinement by computing the percentage error between the numerical solution and the asymptotic value. Roache [20] proposed a GCI method based on Richardson's extrapolation, but this method can be used only in the case of monotonically increasing or decreasing solutions with grid refinement. To overcome this disadvantage, Tanaka and Miyake [21] proposed the least square-based GCI (SLS-GCI) that can be applied in general CFD situations. Their equations describing grid scales are presented below.
where h k is the representative grid size for the three different grids (coarse, medium, and fine), r is the grids ratio, f k is the prediction result of each grid, and is the result difference. The order of convergence (p), extrapolated value ( f c ), and coefficients (α) can be calculated as the solution of Equation (7) as follows: The GCI and error of estimation can be calculated according to the order of convergence.

Bayesian Inference
Bayesian inference is a method to address the inverse problem, where the goal is to back-propagate information about an observation to obtain insight into the uncertainty input parameter. Bayesian inference can calibrate the uncertainty input parameters from the observed data using Bayes' theorem with prior knowledge [22]. where θ are the uncertainty input parameters to be inferred and d is the observation data. p(θ) is the prior distribution of the parameters that reflects the level of information existing on the parameters θ before any measurement is carried out. p(d|θ) is the likelihood (or likelihood function) that represents an estimate of a parameter, and it indicates which parameters are most likely. p(d) is a normalizing factor, known as the evidence or marginal likelihood. p(θ|d) is the posterior distribution of uncertainty for the input parameters.
In the present Bayesian inference, the coefficients of the ZGB cavitation model are considered as uncertainty inputs and are calibrated using experimental observation data by Barre et al. [7]. The overall procedure and framework for Bayesian inference are adopted based on those proposed by Cheung et al. [4].

Point-Collocation Nonintrusive Polynomial Chaos (PC-NIPC)
Bayesian inference includes the process in which the statistical distribution of random inputs is calculated iteratively with a likelihood function and observation data using the specified sampling method. To apply this process in CFD fields, many simulations of the deterministic solver are required to determine the acceptance of the sampling data. When the computational time and cost of each deterministic simulation are considered, the alternative can become the surrogate model or reduced-order statistics analysis. In the present work, PC-NIPC is adopted for the construction of the surrogate model based on the simulation results with sampled inputs [23].
NIPC expansions are functional approximations of a computational model through its spectral representation on a suitably built basis of polynomial functions. The method is to define the joint probability density function (PDF) f X from a random vector with independent components X ∈ R M . The next step is to consider a finite variance computational model as a map Y = M(X), Y ∈ R such that: Then, the polynomial chaos expansion of M(X) is defined as follows: where y α ∈ R are the corresponding coefficients, and A ⊂ N M is the set of selected multi-indices of multivariate polynomials. In Equation (12), ψ α (X) are the multivariate polynomials that are assembled as the tensor product of their univariate counterparts, as shown in Table 2. Table 2. List of classic univariate polynomial families [5].

1(x) 2
Legendre P k (x) In the present study, the PC-NIPC method, which was proposed by Hosder et al. [23], is adopted. For adoption of the PC-NIPC, the minimum sampling number of random input variables is calculated by Equation (13), with the polynomial order (p), the number of random input variables (n), and the oversampling rate (n p ). Hosder et al. [23] recommended that the oversampling rate should be more than 2 for a better approximation with the given polynomial order.
Cheung et al. [4] mentioned that not only is a description of the uncertainty related to the model coefficients but also that a description of the uncertainty regarding the model itself is important. Because the turbulent model coefficients are derived from the assumption of relatively simple flows, such as isotropic turbulence and boundary layer flows, there is a limit to accurately reflecting complex turbulence phenomena in real engineering problems. This model inadequacy becomes a source of uncertainty that propagates through the computational model and leads to errors in the predictions of the model.
In the current work, the model inadequacy of the cavitation model is considered. Assuming d as the true process with no error, the computational model M(X) and the model inadequacy η can be related as follows: Cheung et al. [4] suggested a three-model classification of model inadequacy (no uncertainty, independent Gaussian uncertainty, and correlated Gaussian uncertainty). The correlated Gaussian uncertainty showed the highest posterior plausibility among these model classes [24].

Observation Error
Observation error is one of the unavoidable uncertainties in experimental works. When the parameter uncertainty is calibrated by experimental data corresponding to the observation, modeling the observation error is an important process in Bayesian inference. However, since information on the error bars of the experimental data [7], which are adopted as observations in this simulation, and the data behind these error bars were not provided, the model was assumed to be an independent Gaussian model by Cheung et al. [4]. The true process d, observation error e, and observation data d can be related as follows: where K e is a diagonal matrix with each diagonal component equal to the variance. Observation errors are independent Gaussian random variables with zero mean and standard deviation.

Likelihood Function
All stochastic models are combined into the likelihood function. From Equations (14) and (17), observation data can be expressed by a conditional Gaussian distribution as follows: where σ is the uncertainty parameter including σ e and σ η , and K is the covariance matrix that is the sum of K η + K e . Finally, if all the inferred uncertainty parameters (σ, X) are applied as inputs to find θ, the likelihood function and posterior distribution can be expressed as follows:

Markov Chain Monte Carlo (MCMC)
In practice, posterior distributions do not have a closed-form solution. One widely used option to solve inverse problems relies upon the MCMC model [25]. The basic idea of MCMC simulations is to construct a Markov chain that equals the posterior distribution over the prior support D X . Markov chains can be uniquely defined by their transition probability K x (t+1) x (t) from the step x (t) of the chain at iteration t to the step x (t+1) at the subsequent iteration t + 1. Then, the posterior is equal to the Markov chain if the specified transition probability fulfills the detailed balance condition: One of the algorithms to fulfill this equation is the Metropolis-Hastings (MH) algorithm [26] which is based on proposing and subsequently accepting or rejecting candidate points. At iteration t from the current sample x (t) , one then draws a candidate sample x ( ) from a proposal distribution p x ( ) x (t) . Subsequently, the candidate is accepted with a probability as follows: However, most MCMC algorithms perform poorly when the target distribution shows strong correlations between the parameters. A considerable amount of tuning is necessary to improve the performance of these algorithms. The affine invariant ensemble (AIES) algorithm [27] alleviates this problem by assuming solutions that are invariant with affine transformations of the target distribution. The affine invariance property is achieved by generating proposals according to a so-called stretch move. This refers to proposing a new candidate by: This requires sampling from the distribution p(z) defined by the tuning parameter is then accepted as the new location of the i-th walker, with a probability given as follows:

Geometry and Operating Condition
The two-dimensional geometry is adopted from the experiment of Barre et al. [7]. The height at the venturi throat is 43.7 mm, and the convergence and divergence angles are 4.3 • and 4.0 • , respectively, as shown in Figure 1. The detailed geometric information is presented in Table 3. The conditions of temperature, fluid density, and viscosity of water and vapor are set as shown in Table 4. The velocity at the inlet is set to 10.8 m/s and the pressure at the outlet to 38,450 Pa, shown in Table 5, the same as the experimental condition of Barre et al. [7]. A no-slip boundary condition is applied to the solid wall. The saturated pressure for the prediction of cavitation is given as 2339 Pa.

Geometry and Operating Condition
The two-dimensional geometry is adopted from the experiment of Barre et al. [7]. The height at the venturi throat is 43.7 mm, and the convergence and divergence angles are 4.3° and 4.0°, respectively, as shown in Figure 1. The detailed geometric information is presented in Table 3. The conditions of temperature, fluid density, and viscosity of water and vapor are set as shown in Table 4. The velocity at the inlet is set to 10.8 m/s and the pressure at the outlet to 38,450 Pa, shown in Table 5, the same as the experimental condition of Barre et al. [7]. A no-slip boundary condition is applied to the solid wall. The saturated pressure for the prediction of cavitation is given as 2339 Pa.
From the grid test, which will be explained in Section 3.3, the medium mesh (29,536 nodes) is selected for the present simulation. When it is considered that the mesh size of Barre et al. [7] is 9,861 nodes, the adopted mesh in the present simulation is adequate to simulate the steady sheet cavitation.

Validation of the Deterministic Simulation
Deterministic simulation is validated by the experimental data and numerical simulation result of Barre et al. [7]. Figure 2 shows the comparison of the void fraction and longitudinal velocity (u-velocity) at the stations of X = 5.1 mm, 38.4 mm, and 73.9 mm from the venturi throat. Barre et al. [7] presented two velocity profiles from the experiment; one is the mean velocity (Vmean), and the other is the most probable velocity (Vmp). The most probable velocity (Vmp) is calculated from the probability density function of the measured velocity, which is one of the techniques to predict the velocity in a bubble flow. They showed that mean   From the grid test, which will be explained in Section 3.3, the medium mesh (29,536 nodes) is selected for the present simulation. When it is considered that the mesh size of Barre et al. [7] is 9,861 nodes, the adopted mesh in the present simulation is adequate to simulate the steady sheet cavitation.

Validation of the Deterministic Simulation
Deterministic simulation is validated by the experimental data and numerical simulation result of Barre et al. [7]. Figure 2 shows the comparison of the void fraction and longitudinal velocity (u-velocity) at the stations of X = 5.1 mm, 38.4 mm, and 73.9 mm from the venturi throat. Barre et al. [7] presented two velocity profiles from the experiment; one is the mean velocity (Vmean), and the other is the most probable velocity (Vmp). The most probable velocity (Vmp) is calculated from the probability density function of the measured velocity, which is one of the techniques to predict the velocity in a bubble flow. They showed that mean velocity (Vmean) and the most probable velocity (Vmp) were not the distribution of a vapor bubble in a manner that agrees with the experimental data.
The velocity profiles predicted by the present simulation show a large discrepancy at the station of X = 5.1 mm and 38.4 mm, which is quite similar to the results by Rodio and Congedo [9]. In the simulation by Barre et al. [7], a new extended wall function [28] was adopted to model the turbulent boundary layer, and this model enhanced the prediction of the boundary layer profile relative to the present result. When it is considered that there is also a big difference between two experimental data points, Vmean and Vmp, the prediction of the velocity inside the vapor is still a challenging problem. In the present simulation, one encouraging result is the accurate prediction of the recirculation zone in the wake at the third station (X = 73.9 mm).  The void fraction of a vapor bubble at the wall decreases with an increase in the distance from the venturi throat and is measured as about zero at the station of X = 73.9 mm. The prediction of the void fraction shows good agreement with experimental data at the stations of X = 5.1 mm and 73.9 mm. The thickness and distribution of a vapor bubble show consistent behavior with reference [7] at both positions. However, the simulated void fraction at X = 38.4 mm is similar to that of the numerical simulation by Barre et al. [7]. This overprediction of the void fraction below y~4 mm in two simulation results seems to be due to the limit of multiphase mixture modeling, and this trend is also shown in the simulation results of Rodio and Congedo [9]. In contrast, the simulation result by Barre et al. [7] showed a still large vapor region at the station of X = 73.9 mm; the present simulation is able to predict the thickness and the distribution of a vapor bubble in a manner that agrees with the experimental data.
The velocity profiles predicted by the present simulation show a large discrepancy at the station of X = 5.1 mm and 38.4 mm, which is quite similar to the results by Rodio and Congedo [9]. In the simulation by Barre et al. [7], a new extended wall function [28] was adopted to model the turbulent boundary layer, and this model enhanced the prediction of the boundary layer profile relative to the present result. When it is considered that there is also a big difference between two experimental data points, Vmean and Vmp, the prediction of the velocity inside the vapor is still a challenging problem. In the present simulation, one encouraging result is the accurate prediction of the recirculation zone in the wake at the third station (X = 73.9 mm).
The contour of the void fraction and the longitudinal velocity are presented in Figure 3. As shown in the void fraction contour, as the distance from the venturi throat increases, the thickness of the cavity sheet increases. The length of the sheet, which is determined with the threshold of 0.3 of void fraction, is predicted as 75.0 mm, with a 6.3% difference from the measured value, 80 mm [7]. The velocity contour shows the recirculation zone in the back region of the vapor sheet, which is consistent with the velocity profile at X = 73.9 mm in Figure 2. The contour of the void fraction and the longitudinal velocity are presented in Figure  3. As shown in the void fraction contour, as the distance from the venturi throat increases, the thickness of the cavity sheet increases. The length of the sheet, which is determined with the threshold of 0.3 of void fraction, is predicted as 75.0 mm, with a 6.3% difference from the measured value, 80 mm [7]. The velocity contour shows the recirculation zone in the back region of the vapor sheet, which is consistent with the velocity profile at X = 73.9 mm in Figure 2.

Grid Convergence Index (GCI)
GCI based on the SLS-GCI method [21] is investigated using a coarse, medium, and fine meshes presented in Table 6. The number of nodes in the medium and fine meshes increases by 1.8 and 3.7 times, respectively, compared with the coarse mesh. To evaluate the GCI, the cavitation area, which is calculated with the integration of the cell area with the void fraction below 0.3, is selected as the QoI. The GCI results are summarized in Table  7. The extrapolation values ( ), which mean an estimate of the cavitation area at zero grid spacing, error of estimation ( ), which means the difference between numerical and estimated results, and the uncertainty from the grid ( ) are presented in this table. The cavitation area is extrapolated as 2.53 × 10 −4 at zero grid spacing, and the error and uncertainty between both results are 4.22 × 10 −6 and 1.12 × 10 −5 , respectively. It is judged that the grid test results are reasonable and that the selection of a medium mesh is appropriate in the current work.

Grid Convergence Index (GCI)
GCI based on the SLS-GCI method [21] is investigated using a coarse, medium, and fine meshes presented in Table 6. The number of nodes in the medium and fine meshes increases by 1.8 and 3.7 times, respectively, compared with the coarse mesh. To evaluate the GCI, the cavitation area, which is calculated with the integration of the cell area with the void fraction below 0.3, is selected as the QoI. The GCI results are summarized in Table 7. The extrapolation values ( f c ), which mean an estimate of the cavitation area at zero grid spacing, error of estimation (U e ), which means the difference between numerical and estimated results, and the uncertainty from the grid (U G ) are presented in this table. The cavitation area is extrapolated as 2.53 × 10 −4 at zero grid spacing, and the error and uncertainty between both results are 4.22 × 10 −6 and 1.12 × 10 −5 , respectively. It is judged that the grid test results are reasonable and that the selection of a medium mesh is appropriate in the current work.

Bayesian Inference of ZGB Model Coefficients
In this section, the Bayesian inference of ZGB model coefficients is conducted following the procedure presented in Section 2.3. The first step for Bayesian inference is to determine the random inputs and their prior distribution. In the present work, four coefficients of the ZGB cavitation model are selected as the random inputs, as shown in Table 8. Additionally, two more parameters of model inadequacy and observation error are considered. Because of the lack of knowledge of the model coefficients, it is assumed that these have a uniform distribution. Edeling et al. [29] mentioned that the prior interval of random inputs was selected based on three factors: the spreads of coefficients proposed by the original model, the ranges of these for which the solution was stable, and the avoidance of apparent truncation of the posterior at the edge of the prior range. Since the range of the model coefficients was not presented in the original literature, Zwart et al. [18], the range of each coefficient is initially guessed as ±50% with respect to the original mean initial value. However, since the posterior probability distribution through the calibration procedure is not constrained to the initial prior intervals, the prior range is modified with a wider one or is translated with proper values, as recommended by Edeling et al. [29]. The final ranges of the prior distribution of the four model coefficients are shown in Table 8. The next step is to construct the surrogate model based on PC-NIPC [23] using the sampling data of the random inputs. The calculated sampling number based on PC-NIPC is 30 with a second-order polynomial, two chosen for the oversampling rate, and four random variables.
where K is the covariance matrix of the summation K η model inadequacy and K e observation error [4]. Details regarding the prior distribution are shown in Table 8.
The void fraction profile in the experiment [7] is set to QoI with the assumed error of ±5% and then is adopted as the observation data for Bayesian inference.
To obtain the sample of the posterior distribution, the AIES [27] algorithm is employed with 500 parallel runs of 2000 steps. Each PDF of a parameter was approximated by a kernel density estimation with half of the sample burn-in. After Bayesian inference based on the mentioned procedures, the calibrated model coefficients can be obtained. The posterior distributions of the calibrated coefficients are shown in Figure 4, and the posterior mean values are presented in Table 9.  As shown in Figure 4, the calibrated coefficients are constrained inside the finial prio interval, which means the prior range of the model coefficients is appropriate. The posterior means of ℜ and are decreased from their nominal values, but the posterior means of and are increased after Bayesian inference considering the profile of void fraction as QoI. Since the peak of the void fraction at X = 73.9 mm i predicted to be 75% larger than the original ZGB model coefficients, Bayesian inference adopting the experimental data of the void fraction profile as observation data makes the calibrated model coefficients predict smaller values. The coefficient, , in the ZGB model is related to the production term, and the decrease in this coefficient induce smaller production in the governing Equation (1). The coefficients of and are multiplied in production and destruction terms, respectively. It should be noted that is calibrated as the value 52, larger than the original value of 50, which is a small increase of 4%. In contrast, the coefficient is increased by 30% from 0.01 to 0.013 and then induces an increase in the destruction of cavities.
The void fraction and velocity profiles with the calibrated model coefficients at the station of X = 73.9 mm are presented with the experimental data [7] and the simulation result with the original coefficients in Figure 5. The prediction of the void fraction is 75% better than that of the original prediction, and the velocity profile is still able to correctly  As shown in Figure 4, the calibrated coefficients are constrained inside the finial prior interval, which means the prior range of the model coefficients is appropriate. The posterior means of R B and α nuc are decreased from their nominal values, but the posterior means of F vap and F cond are increased after Bayesian inference considering the profile of void fraction as QoI. Since the peak of the void fraction at X = 73.9 mm is predicted to be 75% larger than the original ZGB model coefficients, Bayesian inference adopting the experimental data of the void fraction profile as observation data makes the calibrated model coefficients predict smaller values. The coefficient, α nuc , in the ZGB model is related to the production term, and the decrease in this coefficient induces smaller production in the governing Equation (1). The coefficients of F vap and F cond are multiplied in production and destruction terms, respectively. It should be noted that F vap is calibrated as the value 52, larger than the original value of 50, which is a small increase of 4%. In contrast, the coefficient F cond is increased by 30% from 0.01 to 0.013 and then induces an increase in the destruction of cavities. The void fraction and velocity profiles with the calibrated model coefficients at the station of X = 73.9 mm are presented with the experimental data [7] and the simulation result with the original coefficients in Figure 5. The prediction of the void fraction is 75% better than that of the original prediction, and the velocity profile is still able to correctly predict the recirculation zone. At the other two stations (X = 5.1 and 38.4 mm), there is a small difference, within 2%, in the void fraction and velocity profiles between both cases. This means that the calibrated model coefficients have little effect on the flow properties of the upstream part. If the velocity profile is adopted as observation data, the model coefficients will be calibrated for a better prediction of the velocity profile. Figure 6 shows the contours of the void fraction with the calibrated model coefficients. It is confirmed that the cavitation area and the sheet length are smaller than the values predicted by the original model.

Uncertainty Quantification with Respect to Operating Conditions
In this section, the uncertainty with respect to the operating conditions in the venturi tube is quantified based on the calibrated ZGB model coefficients. Overall procedures of the corresponding forward problem are (1) the selection of random variables and the determinations of the distributions and their ranges, (2) the definition of QoI, (3) the sampling of random variables depending on the polynomial order and number of inputs, (4) construction of the PC-NIPC model [23], and (5) the calculation of the statistical distribution of QoI and postprocessing.
The considered uncertain input variables are the inlet velocity and the roughness of the wall, which are known as important parameters influencing the flow characteristics of a venturi tube [30,31]. The inlet velocity is modeled as a Gaussian distribution with 10% variation with respect to the given condition, 10.8 m/s. The wall roughness is assumed to be uniform from 0.0001 m to 0.001 m, as shown in Figure 7. In the present work, a wall roughness of zero, an ideally smooth wall, is not considered from the engineering point of view. The QoI is set to the cavitation area for direct comparison with the previous simulation result. The sampling number is determined from 12 to 112, considering the polynomial order from second to sixth orders and the oversampling rate from 2 to 4 for the two input variables and one QoI. The Latin hypercube sampling methodology is adopted to sample random variables, and the PC-NIPC model is constructed for the surrogate model depending on the polynomial order. The leave-one-out (Loo) error, which is an indicator of the proper polynomial order and sampling number, has a minimum of 2.13 × 10 −2 at the fourth polynomial order and an oversampling rate of 4.

Uncertainty Quantification with Respect to Operating Conditions
In this section, the uncertainty with respect to the operating conditions tube is quantified based on the calibrated ZGB model coefficients. Overall the corresponding forward problem are (1) the selection of random vari determinations of the distributions and their ranges, (2) the definition o sampling of random variables depending on the polynomial order and num (4) construction of the PC-NIPC model [23], and (5) the calculation of distribution of QoI and postprocessing.
The considered uncertain input variables are the inlet velocity and the the wall, which are known as important parameters influencing the flow cha a venturi tube [30,31]. The inlet velocity is modeled as a Gaussian distribu variation with respect to the given condition, 10.8 m/s. The wall roughness be uniform from 0.0001 m to 0.001 m, as shown in Figure 7. In the presen roughness of zero, an ideally smooth wall, is not considered from the engi

Uncertainty Quantification with Respect to Operating Conditions
In this section, the uncertainty with respect to the operating conditions in the venturi tube is quantified based on the calibrated ZGB model coefficients. Overall procedures of the corresponding forward problem are (1) the selection of random variables and the determinations of the distributions and their ranges, (2) the definition of QoI, (3) the sampling of random variables depending on the polynomial order and number of inputs, (4) construction of the PC-NIPC model [23], and (5) the calculation of the statistical distribution of QoI and postprocessing.
The considered uncertain input variables are the inlet velocity and the roughness of the wall, which are known as important parameters influencing the flow characteristics of a venturi tube [30,31]. The inlet velocity is modeled as a Gaussian distribution with 10% variation with respect to the given condition, 10.8 m/s. The wall roughness is assumed to be uniform from 0.0001 m to 0.001 m, as shown in Figure 7. In the present work, a wall roughness of zero, an ideally smooth wall, is not considered from the engineering point of view. The QoI is set to the cavitation area for direct comparison with the previous simulation result. The sampling number is determined from 12 to 112, considering the polynomial order from second to sixth orders and the oversampling rate from 2 to 4 for the two input variables and one QoI. The Latin hypercube sampling methodology is adopted to sample random variables, and the PC-NIPC model is constructed for the surrogate model depending on the polynomial order. The leave-one-out (Loo) error, which is an indicator of the proper polynomial order and sampling number, has a minimum of 2.13 × 10 −2 at the fourth polynomial order and an oversampling rate of 4. Figure 8 shows the histogram of the cavitation area, which is the QoI in the present section, according to the variation of operating conditions. The distribution of the cavitation area is similar to an exponential distribution which peaks at areas close to the order of 10 −6 m 2 . The mean of the cavitation area (6.063 × 10 −6 m 2 ) is reduced by 42 times compared with the deterministic solution (2.52 × 10 −4 m 2 ) for the fixed velocity and smooth wall. ergies 2022, 15, 4204 14 o Additionally, the sensitivity of input variables is analyzed using the Sobol index. Sobol indices of the inlet velocity and wall roughness are calculated as 0.59 and 0 respectively. This means that the roughness has much more effect on the cavitation a the QoI, than the inlet velocity.    Additionally, the sensitivity of input variables is analyzed using the Sobol index. The Sobol indices of the inlet velocity and wall roughness are calculated as 0.59 and 0.72, respectively. This means that the roughness has much more effect on the cavitation area, the QoI, than the inlet velocity.   Figure 9 shows the profile of the void fraction and the longitudinal velocities at the two stations of X = 2.5 mm and 5.1 mm. The mean and median values of both are respectively presented in black and orange lines, and the dashed line corresponds to the 95% confidence bound as a result of the uncertainty quantification. In the void fraction profile, there are large 95% confidence bounds starting from 0 to 0.9. However, the velocity profiles show a 95% confidence interval within 10% of the mean value. This means that the wall roughness has a very large effect on the ZGB cavitation model, and the cavitation behavior is affected by the wall roughness. Further investigation of the uncertainty quantification between the cavitation model and wall roughness seems to be necessary. The median area of the void fraction at X = 5.1 mm is predicted to be as small  Additionally, the sensitivity of input variables is analyzed using the Sobol index. The Sobol indices of the inlet velocity and wall roughness are calculated as 0.59 and 0.72, respectively. This means that the roughness has much more effect on the cavitation area, the QoI, than the inlet velocity. Figure 9 shows the profile of the void fraction and the longitudinal velocities at the two stations of X = 2.5 mm and 5.1 mm. The mean and median values of both are respectively presented in black and orange lines, and the dashed line corresponds to the 95% confidence bound as a result of the uncertainty quantification. In the void fraction profile, there are large 95% confidence bounds starting from 0 to 0.9. However, the velocity profiles show a 95% confidence interval within 10% of the mean value. This means that the wall roughness has a very large effect on the ZGB cavitation model, and the cavitation behavior is affected by the wall roughness. Further investigation of the uncertainty quantification between the cavitation model and wall roughness seems to be necessary. The median area of the void fraction at X = 5.1 mm is predicted to be as small as one-third of the mean value, which means that the distribution of the cavitation area resembles the exponential function.

Conclusions
In the present work, Bayesian inference of the cavitation model coefficients is conducted in the venturi tube flow, and the calibrated model coefficients are applied to the venturi tube The contours of the void fraction and longitudinal velocity are presented in Figure 10. It is confirmed that the length of the cavitation sheet is shortened before X = 38.4 mm and that a recirculation zone is shown near X = 38.4 mm.

Conclusions
In the present work, Bayesian inference of the cavitation model coefficients is conducted in the venturi tube flow, and the calibrated model coefficients are applied to the venturi tube flows, which have the same geometry and flow condition as those in Barre et al. [7].
To validate the deterministic solution, the profiles of the void fraction and longitudinal velocity at three stations are compared with experimental data and simulation results by Barre et al. [7]. The void fraction profiles agree well with the reference ones at the front two stations and show better prediction than the simulation

Conclusions
In the present work, Bayesian inference of the cavitation model coefficients is conducted in the venturi tube flow, and the calibrated model coefficients are applied to the venturi tube flows, which have the same geometry and flow condition as those in Barre et al. [7].
To validate the deterministic solution, the profiles of the void fraction and longitudinal velocity at three stations are compared with experimental data and simulation results by Barre et al. [7]. The void fraction profiles agree well with the reference ones at the front two stations and show better prediction than the simulation result by Barre et al. [7]. The prediction of the turbulent boundary layer inside and outside cavitating zone fails, which is consistent with Rodio and Congedo's [9] simulation results, but the recirculation zone is predicted correctly.
GCI analysis by the SLS-CGI method [21] is conducted with respect to the global parameters and cavitation area and shows 2.53 × 10 −4 for the extrapolation value (f c ), 4.22 × 10 −6 for the estimation error (U e ), and 1.12 × 10 −5 for the uncertainty error (U G ). In addition, it is confirmed that the adopted medium mesh with a value approximately equal to the cavitation area of 2.52 × 10 −4 is appropriate.
For Bayesian inference of ZGB model coefficients, four model coefficients are calibrated, and the covariance matrix, which is proposed by Cheung et al. [4], including model inadequacy and observation error, is considered. After initial guesses of the uniform prior intervals of (±50%) of each model coefficient, the range of each is modified with the criteria that the posterior range by Bayesian inference should be constrained inside the prior one. Posterior means of R B and α nuc are decreased from their nominal values, but posterior means of F vap and F cond are increased to decrease the void fraction at the third station at X = 73.8 mm. The velocity profile at all stations shows little change after Bayesian inference. However, there are still large discrepancies in the prediction of velocity and void fraction after observation data are accounted for by Bayesian inference. This means that the perturbation of cavitation model coefficients is not sufficient to remedy the model deficiency. In future work, the other parameters, such as turbulence model inadequacy and discretization error, should be considered simultaneously as uncertain inputs for Bayesian inference [32,33].
Finally, the forward problem with respect to the operating conditions, including the inlet velocity and the wall roughness, is conducted based on the PC-NIPC model with a fourth-order polynomial and an oversampling rate of 4. Therefore, 60 cases are sampled with the Latin hypercube sampling method. When the uncertainty of the wall roughness is considered, in contrast with the previous general simulation with a smooth wall, the cavitation area is decreased 42 times, and the variation of the void fraction is very large, ranging from 0 to 0.9. The discrepancy between the mean and median value of the void fraction is observed. Through the sensitivity analysis with the Sobol index, the wall roughness (0.72) has more effect on the simulation result than the inlet velocity (0.59). This means that the wall roughness in the venturi-type components of the in-service testing should be considered carefully to obtain consistent results with the experiment and quantify the uncertainty of performance in the venturi tube.