X Merle

P Cinnella
email: paola.cinnella@ensam.eu

Robust prediction of dense gas flows under uncertain thermodynamic models

A Bayesian approach is developed to quantify uncertainties associated with the thermodynamic models used for the simulation of dense gas flows, i.e. flows of gases characterized by complex molecules of moderate to high molecular weight, in thermodynamic conditions of the general order of magnitude of the liquid/vapor critical point. The thermodynamic behaviour of dense gases can be modelled through equations of state with various mathematical structures, all involving a set of material-dependent coefficients. For several organic fluids of industrial interest abundant and high-quality thermodynamic data required to specify such coefficients are hardly available, leading to undetermined levels of uncertainty of the equation output. Additionally, the best choice for the kind of equation of state (mathematical form) to be used is not always easy to determine and it is often based on expert opinion. In other terms, equations of state introduce both parametric and model-form uncertainties, which need to be quantified to make reliable predictions of the flow field. In this paper we propose a statistical inference methodology for estimating both kinds of uncertainties simultaneously. Our approach consists of a calibration step and a prediction step. The former allows to infer on the parameters to be input to the equation of state, based on the observation of aerodynamic quantities like pressure measurements at some locations in the dense gas flow. The subsequent prediction step allows to predict unobserved flow configurations based on the inferred posterior distributions of the coefficients. Model-form uncertainties are incorporated in the prediction step by using a Bayesian model averaging (BMA) approach. This consists in constructing an average of the predictions of various competing models weighted by the posterior model probabilities. Bayesian averaging also provides a useful tool for making robust predictions from a set of alternative calibration scenarios (Bayesian model-scenario averaging or BMSA). The proposed methodology is assessed for a class of dense gas flows, namely transonic flows around an isolated airfoil, at various free-stream thermodynamic conditions in the dense-gas region.

Introduction

Flows of dense gases, i.e. flows of organic fluids of moderate to high molecular weight working close to saturation conditions, are encountered in several engineering problems, one of the most attractive applications being represented by energy conversion cycles, like heatpumps, refrigeration and, most of all, Organic Rankine Cycles [START_REF] Monaco | Supersonic flows of dense gases in cascade configurations[END_REF][START_REF] Brown | Application of bethe-zel'dovich-thompson fluids in organic rankine cycle engines[END_REF][START_REF] Horen | Numerical simulation of real-gas flow in a supersonic turbine nozzle ring[END_REF][START_REF] Kirillov | Analysis of modern natural gas liquefaction technologies[END_REF][START_REF] Zamfirescu | Performance investigation of high-temperature heat pumps with various bzt working fluids[END_REF].

Dense gas flow simulations can be extremely sensitive to the model used to describe the fluid thermodynamic behavior and its closure coefficients [START_REF] Cinnella | Sensivity analysis of dense gas flow simulations to thermodynamic uncertainties[END_REF], i.e. to the equations of state (EOS). Specifically, EOS give raise to two kinds of uncertainties: the first one concerns choosing a suitable mathematical form among the many available (e.g., cubic EOS [7][START_REF] Soave | Equilibrium constants from a modified redlich-kwong equation of state[END_REF][START_REF] Stryjek | An improved peng-robinson equation of state for pure compounds and mixtures[END_REF], virial EOS with a more or less large number of expansion terms [START_REF] Martin | Development of an equation of state for gases[END_REF][START_REF] Benedict | An empirical equation for thermodynamic properties of light hydrocarbons and their mixtures: i. methane, ethane, propane, and nbutane[END_REF], reference EOS based on power-law expansions of the Helmoltz free energy [START_REF] Span | Equations of state for technical applications. i. simultaneously optimized functional forms for non polar and polar fluids[END_REF]); on the other hand, the material-dependent coefficients associated to the EOS are often imperfectly known, especially for complex organic fluids for which abundant high quality data are less readily available than for widely employed light gases like hydrogen, nitrogen, carbon dioxide, etc [START_REF] Lemmon | Short fundamental equations of state for 20 industrial fluids[END_REF]. Although parametric uncertainty may be critical for the accurate prediction of dense gas flows, in [START_REF] Cinnella | Sensivity analysis of dense gas flow simulations to thermodynamic uncertainties[END_REF] it was shown that for some complex gases the model-form uncertainty can be even overwhelming with respect to the former one.

Several studies have addressed the problem of calibrating EOS from thermodynamic data available for the fluid(s) of interest. This can be done either in a deterministic way (e.g. [START_REF] Lemmon | Short fundamental equations of state for 20 industrial fluids[END_REF] and references cited therein) or stochastically, e.g. by using statistical inference methodologies [START_REF] Hajipour | Uncertainty analysis applied to thermodynamic models and process design -i. pure components[END_REF][START_REF] Robinson | Fundamental issues in the representation and propagation of uncertain equation of state information in shock hydrodynamics[END_REF]. In a previous work [START_REF] Merle | Bayesian quantification of thermodynamic uncertainties in dense gas flows[END_REF], the present authors used aerodynamic data instead of thermodynamic ones for calibrating the material-dependent coefficients of equations of state (EOS) used to model the thermodynamic behavior of the working fluid. For that purpose, a Bayesian statistical procedure was used to infer on the posterior probability distributions of the closure coefficients associated with three well-known thermodynamic models, given data on the wall pressure distribution for a dense gas flow past an airfoil. The statistical model used for the calibration accounted both for uncertainties in the observed data and in the model form. The latter expresses the fact that, due to the simplifying assumptions intrinsic to any mathematical model of a physical system, this can never predict exactly the observed values, even assuming that the best possible coefficients are available (see [START_REF] Brynjarsdottir | Learning about physical parameters: the importance of model discrepancy[END_REF] for a thorough discussion on the role of model-form uncertainties). The model-form uncertainty was represented as a Gaussian random vector with a given correlation structure, whose coefficients (called the hyperparameters) were calibrated alongside the physical model parameters.

Accounting for model-form uncertainty in the statistical calibration is useful to temperate overfitting problems (see [START_REF] Merle | Bayesian quantification of thermodynamic uncertainties in dense gas flows[END_REF][START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF]). Nevertheless, such problems are not completely avoided and the closure parameter posterior distributions have limited validity when used to predict flow configurations far away the one on which they where calibrated. In other terms, they have no universal validity [START_REF] Brynjarsdottir | Learning about physical parameters: the importance of model discrepancy[END_REF]. Additionally, the model-form uncertainty term is specific to the kind of data used for the calibration and cannot be used, e.g., to predict a different quantity of interest (QoI) for a new flow.

A coherent framework for making predictions in situations where multiple competing models are available is represented by multi-model approaches, used in a plethora of applications including oil price predictions, meteorology, ground-water modeling, aerodynamics, and aeroelasticity [START_REF] Draper | Assessment and propagation of model uncertainty[END_REF][START_REF] Poroseva | Improving the predictive capability of turbulence models using evidence theory[END_REF][START_REF] Duan | Multi-model ensemble hydrologic prediction using bayesian model averaging[END_REF][START_REF] Tebaldi | The use of the multi-model ensemble in probabilistic climate projections[END_REF][START_REF] Diomede | The use of the multi-model ensemble in probabilistic climate projections[END_REF][START_REF] Rojas | Application of a multimodel approach to account for conceptual model and scenario uncertainties in groundwater modelling[END_REF][START_REF] Nitschke | Model-form and predictive uncertainty quantification in linear aeroelasticity[END_REF][START_REF] Raftery | Using bayesian model averaging to calibrate forecast ensembles[END_REF]. Bayesian model averaging (BMA) [START_REF] Draper | Assessment and propagation of model uncertainty[END_REF][START_REF] Hoeting | Bayesian model averaging: a tutorial[END_REF] is among the most widely used multi-model approaches, where posterior average predictions are inferred by weighing individual forecasts from competing models based on their relative skill, with predictions from better performing models receiving higher weights than those of worse performing models. BMA avoids having to choose a model over the others and provides instead a measure of the model-form uncertainty based on the level of agreement among the competing models considered in the average. Due to the high computational cost when combined with complex computer models, BMA has found application in computational fluid dynamics problems, governed by complex nonlinear equations, only recently. In [START_REF] Edeling | Predictive rans simulations via bayesian modelscenario averaging[END_REF] a BMA approach was used to investigate model-form uncertainty due to the existence of several competing closure models for the Reynolds-Averaged Navier-Stokes (RANS) equations, often used to describe turbulent flows. More specifically, an extended formulation of BMA was adopted (first suggested in [START_REF] Draper | Assessment and propagation of model uncertainty[END_REF], see also [START_REF] Rojas | Application of a multimodel approach to account for conceptual model and scenario uncertainties in groundwater modelling[END_REF] for an application to ground-water modeling), which accounts also for the uncertainty about the validity of the calibrated model parameters when applied to a new prediction scenario. Such an extension, termed Bayesian Model-Scenario Averaging (BMSA) in [START_REF] Edeling | Predictive rans simulations via bayesian modelscenario averaging[END_REF], amounts to considering each realization of a model supplemented by coefficients adjusted on different calibration scenarios as a component of the mixture, weighted through a suitable 'scenario' probability. In this framework, the posterior predictive distribution of a QoI for a model applied to a new scenario is the weighted average of the predictions of a model using different sets of coefficients. An attractive feature of BMSA is that, if new sets of coefficients become available, these can easily be added to the BMSA mixture, which is then expanded to a larger set of scenarios. Differently from BMA, where the model weights are inferred during the calibration process, in BMSA the scenario weights have to be chosen a priori, according to expert judgement or Bayesian criteria. Eventually, if data are available for the new scenario, these can be used to infer on the scenario weights. This is however not the case in general, so that a key point is to find a suitable a priori criterion for choosing the scenario weights. An empirical criterion, based on the level of agreement among the competing models calibrated on the same scenario, was proposed in [START_REF] Edeling | Predictive rans simulations via bayesian modelscenario averaging[END_REF] and successfully applied to the prediction of a variety of turbulent perfect gas flows in [START_REF] Edeling | Predictive rans simulations via bayesian modelscenario averaging[END_REF][START_REF] Edeling | Bayesian predictions of reynoldsaveraged navier-stokes uncertainties using maximum a posteriori estimates[END_REF].

Accounting for calibration scenario uncertainty is an essential feature for dense gas flows, which are very sensitive to the thermodynamic operating conditions (see, e.g., [START_REF] Cinnella | Bayesian quantification of thermodynamic uncertainties in dense gas flows[END_REF]). A standard approach to calibrate a model over data sets coming from different scenarios consists in inferring from all data sets simultaneously. For instance, Cheung et al. [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF]. inferred on the coefficients of a turbulence model using velocity and skin friction data for flat plate boundary layers subject to three different external pressure gradients. In such an approach, the resulting coefficients represent a compromise allowing the model to fit reasonably well all of the available data. The calibration procedure can than be applied to different models, and a final BMA can be used to predict a new case while accounting for model-form uncertainty. The interest of such an approach is that it avoids choosing scenario weights. The drawback is that the likelihood function involved in the inference of such a large data set is a multivariate random function, with a highdimensional correlation matrix that may be hard to build and invert. Additionally, including new data sets implies restarting the calibration from scratch. For this reason, in this paper we propose an alternative approach, termed Bayesian Scenario Averaging (BSA). This can be seen as a special case of BMSA, where a single model is calibrated against several competing scenarios and applied to a new one. A suitable a priori criterion is proposed to weight the scenarios. BSA is assessed against the standard approach in the case of dense gas flow predictions.

The main goal of the present work is to develop a robust tool for prediction of dense gas flows with quantified thermodynamic uncertainties, based on Bayesian inference and model averaging. Specifically, we first use Bayesian calibration for updating the coefficients of two different EOS, namely the cubic EOS of Peng-Robinson-Strjyek-Vera [START_REF] Stryjek | An improved peng-robinson equation of state for pure compounds and mixtures[END_REF] (PRSV) and the 5th-order virial Martin-Hou EOS [START_REF] Martin | Development of an equation of state for gases[END_REF] (MAH). For that purpose, similarly to [START_REF] Merle | Bayesian quantification of thermodynamic uncertainties in dense gas flows[END_REF], we consider pressure data for dense gas flows past an airfoil [START_REF] Rusak | Transonic flow of dense gases around an airfoil with a parabolic nose[END_REF][START_REF] Wang | Numerical studies of transonic BZT gas flows around thin airfoils[END_REF][START_REF] Rusak | Low-drag airfoils for transonic flow of dense gases[END_REF] characterized by different free-stream thermodynamic conditions. For all the flows, the airfoil geometry and the free-stream Mach number are the same, while the free-stream pressure and temperature may vary, leading to a different thermodynamic behaviour. In this context, we term a scenario a given choice of the operating thermodynamic conditions. The EOS are calibrated versus different scenarios individually or simultaneously. In the first case (called single-point calibration), posteriors of the coefficients resulting from different calibrations are propagated through a new scenario by means of BSA. The advantages of the BSA approach and the role of the weighting criterion are discussed for two prediction scenarios (not used for calibration) for which validation data is available. Finally, information from different scenarios and models is mixed together using BMSA, demonstrating the importance of accounting for both parametric and model-form uncertainty.

The paper is organized as follows: in Section 2 we present the physical problem, the governing equations and the thermodynamic models (EOS), and we briefly describe the numerical solver used for the computations. In Section 3 we describe the data used for the calibrations and the choice of the calibration and prediction scenarios. Section 4 illustrates the Bayesian calibration methodology and the BMSA framework for robust prediction. Numerical applications of the calibration and prediction approach to transonic dense gas flows are reported in Sections 5 and 6. Finally, Section 7 is devoted to conclusions and perspectives.

Description of the problem

The objective of this study is twofold: (i) first we compare two calibration strategies to take into account data coming from multiple flow scenarios, (ii) we assess the capability of BMSA to provide robust predictions of dense gas flows. Such flows are of special interest for energy conversion machines which involve an expander (often a turbine) and/ or a compressor, characterized by the presence of bladed disks, through which the dense gas flows. Since blade sections can be roughly seen as airfoils, hereafter we investigate the feasibility of our calibration methodology for a simplified configuration, roughly representative of a blade section, i.e. an isolated airfoil. Such kind of configuration was used in the past to investigate qualitatively dense gas effects in turbomachinery [START_REF] Cramer | Transonic flows of bethe-zel'dovich-thompson fluids[END_REF][START_REF] Cinnella | Optimal airfoil shapes for viscous transonic flows of bethe-zel'dovich-thompson fluids[END_REF]. Moreover, for this simplified problem sensitivity studies of the impact of thermodynamic uncertainties on flow simulation results where carried out in [START_REF] Cinnella | Sensivity analysis of dense gas flow simulations to thermodynamic uncertainties[END_REF].

Precisely, the case selected for this study is the steady transonic flow of a dense gas over a thin wing section, namely, a NACA0012 airfoil. The working fluid is a siloxane (silicon oil), known with the commercial name of D5 (chemical formula ((CH 3 ) 2 SiO) 5 ).

Governing equations

For simplicity, the flow around the airfoil is assumed to be inviscid, which represents a good approximation as long as the QoI is represented by the wall pressure distribution and the flow boundary layer remains attached. This amounts to consider that the flow is governed by the compressible Euler equations, written hereafter in integral form for a control volume Ω with boundary ∂Ω:

+ = w f n d dt d dS • 0 (1)
In Eq. ( 1), w is the conservative variable vector, where

= w v E ( , ,
) T n is the outer normal to ∂Ω, and f, is the flux density:

= f v vv v I H ( , ¯, , ) T
where ρ is the density, v is the velocity vector, E the specific total energy, = + H E p/ the specific total enthalpy, p is the pressure and I ¯is the unit tensor. The preceding equations are completed by a thermal equation of state:

= w w p p T ( ( ), ( )) (2)
with T the total pressure, and by a caloric equation of state for the specific internal energy e, which must satisfy the compatibility relation:

= = + w w e e T e c T dT T p T p d ( ( ), ( )) ( ) T T v 0 , 2 0 0 (3)
In Eq. ( 3), c v, ∞ is the ideal gas specific heat at constant volume, quantities with a prime superscript are auxiliary integration variables, and subscript 0 indicates a reference state. The caloric equation of state is completely determined once a variation law for c v, ∞ has been specified.

The problem setup is completed by specifying the free-stream conditions, i.e. flow Mach number and angle of attack and thermodynamic conditions at the far field. These are treated in the following as deterministic and equal to = M 0.95 and = AoA 0 0 for all cases. On the other hand, the free-stream thermodynamic conditions, given in terms of reduced pressure = p p p / r c , and density = / , r c , are used to specify the different flow configurations, i.e. the scenarios.

For dense gases, advanced EOS are needed to account for their complex thermodynamic behaviour, which may be described by means of the fundamental derivative of gas dynamics [START_REF] Thompson | A fundamental derivative in gas dynamics[END_REF] = + a a : 1

s (4)
where ρ is the fluid density, a the sound speed, and s the entropy. Γ can be interpreted as a measure of the rate of change of the sound speed with density in isentropic perturbations. If Γ < 1, the flow exhibits an uncommon sound speed variation in isentropic perturbations: a grows in isentropic expansions and decreases in isentropic compressions, the opposite of what happens in "common" fluids. For instance, in perfect gases, Γ is equal to + ( 1)/2, where the specific heats ratio γ is always greater than 1 for thermodynamic stability reasons; therefore, Γ > 1 as well. Finally for some heavy polyatomic fluids, referred to as the Bethe-Zel'dovich-Thompson (BZT) fluids, Γ may take negative values in a small thermodynamic region above the liquid/vapor saturation curve, leading to non classical behaviors in the transonic and supersonic regime, such as expansion shock waves, mixed shock/fan waves, and splitting shocks (see, e.g. [START_REF] Cinnella | Optimal airfoil shapes for viscous transonic flows of bethe-zel'dovich-thompson fluids[END_REF] and references cited therein). The thermodynamic region in the vapour phase where Γ < 1 is often termed the dense gas region, while the vapor region where Γ < 0 is called the inversion zone. The presence of regions of low or negative values of the fundamental derivative within a flow field may change dramatically the behaviour compared to that of a classical gas. This is why the size and location of the inversion zone has a deep impact on the resulting flow field. The inversion zone has been found to be extremely sensitive to the equation of state in use and the associated input parameters (see e.g. [START_REF] Cinnella | Sensivity analysis of dense gas flow simulations to thermodynamic uncertainties[END_REF]).

In this work, we consider two alternative thermodynamic models for the thermal equation of state, namely, the Peng-Robinson-Stryjek-Vera (PRSV) cubic equation of state [START_REF] Stryjek | An improved peng-robinson equation of state for pure compounds and mixtures[END_REF], and the multi-parameter Martin-Hou (MAH) equation [START_REF] Martin | Development of an equation of state for gases[END_REF], based on a five-term virial expansion. Both models have been often used in the dense gas literature, see e.g. [START_REF] Cinnella | Sensivity analysis of dense gas flow simulations to thermodynamic uncertainties[END_REF][START_REF] Merle | Bayesian quantification of thermodynamic uncertainties in dense gas flows[END_REF].

In addition to the preceding EOS, we also consider for calibration and validation purposes (as discussed later) a more complex and accurate thermodynamic model, namely, the multiparameter technical equation of state introduced by Span and Wagner (SW) [START_REF] Span | Equations of state for technical applications. i. simultaneously optimized functional forms for non polar and polar fluids[END_REF] with coefficients adjusted to D5 in [START_REF] Colonna | Multiparameter equations of state for selected siloxanes[END_REF]. This is a complex technical equation of state, which is considered as the most accurate and complete thermodynamic model presently available for this fluid, provided that a sufficiently large set of high-quality experimental data is available for a reliable fitting of the model coefficients, i.e. for reproducing experimental data with very high accuracy. One drawback is that this EOS is much more costly than the preceding ones. In practice, the equation coefficients are obtained by means of a regression procedure: the experimental data used for calibration are weighted in such a way that high-quality data contribute with higher weights. Experimental data with uncertainties larger than a required threshold are discarded, or, when no or insufficient data are available for a certain region, they are weighted according to their level of uncertainty. Note that, due to the limited and highly uncertain data available for the dense gas of interest, the SW model is in practice also affected by significant uncertainties, of the order of about 10%. In the following calibrations, these are treated as "observational" uncertainties.

The main features of the three models are briefly recalled hereafter. We refer to [START_REF] Cinnella | Sensivity analysis of dense gas flow simulations to thermodynamic uncertainties[END_REF], [START_REF] Merle | Bayesian quantification of thermodynamic uncertainties in dense gas flows[END_REF] for a discussion of their sensitivity to uncertain thermodynamic input parameters.

Thermodynamic models PRSV model

The cubic Peng and Robinson equation of state with modifications suggested by Stryjek and Vera [START_REF] Stryjek | An improved peng-robinson equation of state for pure compounds and mixtures[END_REF] can be written, in reduced form (i.e. with thermodynamic quantities normalised with respect to their values at the liquid/vapour thermodynamic critical point), as:

= + p T Z v b a T v b v b / ( ) 2 r r c r r r r r r r 2 2 (5)
where the subscript r denotes a reduced quantity,

= (•) (•)/(•) ,
The critical compressibility factor is unequivocally determined by imposing that p r equals 1 at the critical point, i.e. for = v 1, ). Thus, the only free parameter left in the thermal equation is the acentric factor ω.

The PRSV equation is supplemented with a model for the ideal gas contribution to the specific heat at constant volume, represented here by a power law of the form:

= c T c T T T ( ) ( ) v v c c n , , (9)
where the c subscript denotes critical-point values, and the exponent n and the ideal-gas-limit isocoric specific heat at the critical temperature c v, ∞ (T c ) are material-dependent constants. Finally, the non-dimensional PRSV model [START_REF] Zamfirescu | Performance investigation of high-temperature heat pumps with various bzt working fluids[END_REF][START_REF] Stryjek | An improved peng-robinson equation of state for pure compounds and mixtures[END_REF] depends on the uncertain parameters ω, n, and c v, ∞ (T c ).

MAH model

The comprehensive thermal equation of state of Martin and Hou [START_REF] Martin | Development of an equation of state for gases[END_REF] reads:

= + = p T Z v b f T v b ( ) ( ) ( ) r r c r r i r i r r r i 2 5 , (10)
= + + f T A B T C kT ( ) exp( ) r i r r i r i r r i r , , , ,
with = k 5.475. The gas-dependent coefficients A r, i , B r, i , C r, i can be expressed in terms of the critical temperature and pressure, the critical compressibility factor, the Boyle temperature (which may be expressed as a function of the critical temperature) and one point on the vapour pressure curve. The MAH equation of state is supplemented again by Eq. ( 9) to compute the ideal gas contribution to the specific heat at constant volume. Globally, the MAH thermodynamic model in reduced form, Eqs. ( 10) and ( 9), requires the knowledge of six material-dependent parameters, namely, the critical pressure p c , the critical temperature T c , the critical compressibility factor Z c , the normal boiling temperature T b , the exponent n and the reduced ideal-gas constantvolume specific heat at the critical temperature c v, ∞ (T c ). As a consequence, we choose to neglect them for the following of the study. Nominal values of the different input parameters corresponding to D5 are given in Table 1.

SW model

The last thermodynamic model considered in this study is a 12parameter technical equation of state based on the functional form for non polar fluids proposed by Span and Wagner ([12]). This is written as an expression for the reduced Helmholtz free energy Φ (i.e. normalised with RT c ), sum of an ideal-gas part, Φ 0 , function of the ideal-gas isobaric heat capacity c p, ∞ , and by a residual term Φ r that takes into account real-gas corrections:

= + + + + + + + + + + + n n n n n n n n n n n n ( , ) exp( ) exp( ) exp( ) exp( ) exp( ) exp( ) exp( ) r 1 0.25 2
where n 1 , ... , n 12 are substance-specific coefficients, = / c is the reduced density and = T T / c is the inverse of the reduced temperature. Material-dependent coefficients for D5 have been taken from Ref [START_REF] Colonna | Multiparameter equations of state for selected siloxanes[END_REF]., to which we refer for more details. The thermal and caloric EOS are derived from the following thermodynamic relations:

= + = + p RT e RT 1 ; r r 0
For the calculation of caloric properties, the SW EOS (Eq. ( 11)) is supplemented by the ideal gas contribution to the specific heat at constant pressure, which is now approximated here as a polynomial function of the temperature:

= + = + + + c T R c T R c c T c T c T ( ) ( ) 1 p r v r , , 1 2 3 2 4 3 (12)
where the polynomial coefficients c i depend again on the substance under consideration and are given for D5 in [START_REF] Colonna | Multiparameter equations of state for selected siloxanes[END_REF]. This model is considered as more accurate than the preceding ones. For this reason, it is used in the following to generate reference solutions used to calibrate the simpler PRSV and MAH models.

Numerical solver

Numerical solutions for the dense gas flow of interest are found by means of an in-house dense gas flow solver. The governing equations are discretized using a cell-centered finite volume scheme for structured multi-block meshes of nominal third-order accuracy, which allows computing flows governed by an arbitrary equation of state. We refer to [START_REF] Cinnella | Numerical solver for dense gas flows[END_REF] and the references therein for more details about the numerical solver and its validation.

Due to high computational cost, the dense gas solver is used to generate an inexpensive surrogate model, i.e. an analytical function providing an approximation of the code output for a given choice of the input parameters. The surrogate model used in this work is a deterministic piecewise polynomial interpolation. A more detailed description of the surrogate model, along with a discussion about the approximation error can be found in [START_REF] Merle | Bayesian quantification of thermodynamic uncertainties in dense gas flows[END_REF].

Calibration data

Pseudo-experiment

Despite many past efforts for carrying out experiments for dense gas flows [START_REF] Borisov | Rarefaction shock waves near the critic liquid-vapour point[END_REF][START_REF] Fergason | Theory for producing a single-phase rarefaction shock-wave in a shock tube[END_REF][START_REF] Colonna | Design of the dense gas flexible asymmetric shock tube[END_REF], no detailed experimental characterization of dense gas flows is available yet. Very recently, preliminary results have been presented for nozzle flows of the light siloxane MDM [START_REF] Spinelli | Preliminary characterization of an expanding flow of siloxane vapor MDM[END_REF]. These include flow visualizations and a single pressure measurement at the geometrical nozzle throat. Other dense-gas facilities are under development, e.g. [START_REF] Reinker | Thermodynamics and fluid mechanics of a closed blade cascade wind tunnel for organic vapors[END_REF], so that a more complete body of experimental knowledge will come into availability in the near future. In this work, to provide a proof of concept of our predictive methodology, we generate synthetic "experimental" data by running the dense gas solver with the most complex and accurate thermodynamic model, namely, the SW model.

This numerical experiments consists of transonic steady flows of siloxane D5 past an NACA0012 airfoil with = M 0.95, = AoA 0 0 and various choices of the free-stream thermodynamic conditions (corresponding to different scenarios), which are discussed later. Values of the pressure coefficient:

= C p p U p 1 2 2 (13)
(where U ∞ is the free-stream velocity) are collected at 17 locations distributed along the airfoil upper wall (simulated pressure taps) highlighted by red symbols in Fig. 1(a)). A more complete description of the pseudo-experimental setup can be found in [START_REF] Merle | Bayesian quantification of thermodynamic uncertainties in dense gas flows[END_REF]. In Fig. 1(a) we report also typical isoline of the pressure coefficient around the airfoil, showing that the flow is characterized by a shock wave on the rear part of the airfoil for the chosen Mach number, as well as a typical

Table 1

Nominal values of thermodynamic properties and parameters for D5. The collected synthetic pressure data are perturbed by adding a Gaussian noise with zero mean and standard deviation equal to 10% of the nominal value to simulate the effect of experimental errors. This kind of observational error model is often used in the literature [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF]. The 10% value for the standard deviation was chosen based on preliminary investigations of the sensitivity of the Span Wagner thermodynamic model to uncertainties in the caloric quantities [START_REF] Cinnella | Sensivity analysis of dense gas flow simulations to thermodynamic uncertainties[END_REF]).

Choice of the calibration and prediction scenarios

Pseudo-experimental data are generated for three choices of the free-stream thermodynamic conditions close to the liquid/vapour saturation curve, where dense gas effects play a crucial role in the flow physics. Their locations in the Clapeyron p v diagram are depicted in Figure 2 and the corresponding values of the reduced thermodynamic conditions (numbered from 1 to 3) are given in Table 2. Furthermore, pseudo-experimental data were also generated for two additional operating conditions, not used in the calibration process, corresponding to conditions 4 and 5 in Table 2 and reported with green symbols in Fig. 2. Such conditions are referred to as the prediction scenarios, and the corresponding data are used to validate the predictive model. In this figure we also report the saturation curve computed according to the reference SW model as well as the saturation curves predicted by the PRSV and MAH models with parameters set to their nominal values. In all cases, the operating conditions are sufficiently far from the saturation curve to avoid the appearance of two-phase flow conditions in the calculation. We checked that this condition was satisfied also for the perturbed values of the parameters considered in the calibration and prediction procedures.

Bayesian methodology

Bayesian calibration

This section describes the statistical procedure used to calibrate the closure parameters associated to the EOS described in Section 2.

We observe that a QoI y can be computed as an output of the model M given a set of parameters n : To that aim, consider first a generic physical problem of the form:

= y M ( ) (14)
Note that y may also depend on additional parameters (explanatory variables) that are considered as known and do not need to be calibrated (for instance the free-stream conditions or the airfoil geometry). Such fixed parameters are not explicitly represented as an argument for M.

In the Bayesian framework, the unknown parameter vector θ is treated as a random vector, characterised by means of its joint probability density function (pdf), noted as f. Due to the uncertainty on θ, y is a dependent random vector.

The main goal of Bayesian calibration is to achieve new knowledge about θ by constructing an improved representation of its pdf, starting from some prior knowledge on θ and some observed data.

For this purpose, assuming that new amount of information is available and represented by a random vector of observed data D , N the Bayes rule simply states that:

= = = = f D D f D D f D D f ( | ) ( | ) ( ) ( ) (15)
where:

• f(θ) is the prior pdf and represents the initial belief about θ;

• = f D D ( ¯| )
is the likelihood, and corresponds to the probability to observe D ¯, a realisation of the random variable D, if θ is known exactly;

• = f D D ( | ¯)
is the posterior pdf and represents the updated knowledge of θ given information about D;

• = f D D ( ¯)
is the evidence, which represents the probability to observe D ¯for all the values of θ, and acts as a normalisation constant so that Eq. ( 15) is shortly written:

= = f D D f D D f ( | ) ( | ) ( ). (16)
Calibration is based on comparisons between modelled and observed data. Thus, it is necessary to extract from all possible model outputs quantities corresponding to the observed data. Note that the calibrated model can subsequently be used to predict other QoI, i.e. the output vector y includes, but is not necessarily equal to D. In the following of this section we assume for simplicity that = y D. To make this point appear in the former notation, we rewrite Eq. ( 16):

= = f D D y f D D y f ( | , ) ( | , ) ( ). (17)
even if this notation is redundant. For our concerns, the uncertain parameter vector θ is the set of closure parameters associated to a given EOS, y is a vector of outputs from the dense gas solver M (in our case, the computed pressures at the 17 locations for which data are available), and D is a random vector of observed data (i.e. the reference pressure values).

From Eq. ( 17), it appears that the distribution of the prior and the likelihood have a strong influence on the posterior and have to be chosen carefully. A number of priors have been proposed and used in the literature. Ref. [START_REF] Arnst | Predictive rans simulations via bayesian model-scenario averaging[END_REF] recommends using the improper uniform pdf if no information is available on the coefficients in advance of observing the data, and a maximum entropy prior [START_REF] Soize | Construction of probability distributions in high dimension using the maximum entropy principle: applications to stochastic processes, random fields and random matrices[END_REF] when information concerning their mean, covariance or other generalized moments is available. In the context of the present paper we only expect the uncertain parameters to be positive and finite, but we do not have reasons to enforce any preferential parameter range. This is why we chose for the priors non-informative uniform distributions:

f a b ( ) ([ , ]) (18)
where the interval [a, b] is taken large enough to ensure a good exploration of the parameter space, while avoiding to include nonphysical values. The likelihood function is a statistical model describing the probability of observing the data for a given choice of the model parameters. It may include information both about the observational error on the data and the model error, including the effect of surrogate modelling if needed. In the following, we use a likelihood function resulting from a multiplicative/additive error model, similar to that used in [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF] for calibrating the closure coefficients of turbulence model from measured velocity profiles in a turbulent boundary layer. We refer to [START_REF] Kennedy | Bayesian calibration of computer models[END_REF], [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF] and the references cited therein for more details about possible choices for the construction of likelihood functions.

In this work, the model includes an additive error modelling observational noise and a multiplicative term representative of the model inadequacy. First, the data D at a given location x i along the airfoil wall are modelled as:

= + D x D x e ( ) ^( ) i i i ( 19
)
with

• e i , the experimental noise at location x i ;

• D x ^( ), i the true (unobserved) pressure coefficient value at x i .

The observational errors,

= … e i N ,
1, , i ( = N 17 being the number of data) are taken here as independent and normally distributed with zero mean and a standard deviation equal to 10% of the observed value. On the other hand, the true process D ^is in general not exactly captured by the model ( 14) -even with the best possible choice of the parameters θ -because of inadequacies intrinsic to the simplifying assumption used to mathematically modelling the physical phenomenon under study. In this work, we assume the true process to be equal to the model output y (x i , θ), multiplied by an error coefficient η i :

= D x y x ^( ) ( , ) i i i ( 20
)
which takes into account the discrepancy between the simulation and the actual system. Note that the use of corrective multiplicative coefficients to account for model predictive deficiencies (safety factors) is common engineering practice (see, e.g. [START_REF] Ullman | The mechanical design process[END_REF]). The vector of model-inadequacy terms = … ( , , ) N 1 is assumed to be well represented by a correlated Gaussian model of the form:

K (1, ) M with = K x x i j N ( ) exp ( ) 2 , 1 , M ij i j 2 2 2 (21)
where x i and x j are two distinct observation abscissas separated by the length scale α. The coefficients σ and α are supplementary parameters intrinsic to the statistical model (hyperparameters) and are calibrated along with the parameters of the physical model. Note that σ represents the magnitude of the model inadequacy and thus can be taken as an indicator of the accuracy of a given model. On the other hand, the model-inadequacy term is expected to reduce the risk of over-fitting the model parameters, by introducing additional degrees of freedom (see, e.g. [START_REF] Brynjarsdottir | Learning about physical parameters: the importance of model discrepancy[END_REF] for a discussion about the importance of model-inadequacy terms in Bayesian calibration). The correlation kernel ( 21) is an adaptation of the Matérn model suggested by Kennedy and O'Hagan [START_REF] Kennedy | Bayesian calibration of computer models[END_REF] and it is an expression of the belief that the code has some degree of smoothness, so that errors at closeby observation points (within a length α) are correlated.

With the preceding assumptions, the product

= … y ( ) i i i N
1, , happens to follow a normal distribution with mean y and a covariance matrix defined by:

= K y y K i j N ( ) ( ) , 1 , M ij i j M ij (22)
Finally, because

= … e ( ) i i N 1, ,
is also a random vector, the likelihood can be written under the form:

= f D y K D y K D y ( | , ) 1 (2 ) | | exp 1 2 ( ) ( ) N T 1 (23)
where

= + K K K , e M
with K e a diagonal matrix associated to the observational error vector.

An interesting point of the Bayesian framework is the possibility to account for various calibration scenarios in a one shot procedure by considering a meta vector of data made of several D ¯j . In such a configuration, the covariance matrix K M must be carefully modified to consider the correlations, not only between the pressure taps within a particular scenario, but also between scenarios. In this work, the covariance matrix accounts for correlations between pressure predictions for different scenarios. As such, considering J different scenarios, K M is now written

= + = + = + × K x x p p v v i j N k l J s k N i t l N j s t J N ( ) exp ( ) 2 exp ( ) ( ) 2 , with 1 , 1 , and ( 1) ( 1) 1
,

M st i j k l k l 2 2 2 2 2 2 ( 24
)
where p and v are the pressure and the specific volume and define a point of operating conditions in Clapeyron's diagram; β is a new hyper parameter corresponding to a correlation length in Clapeyron's diagram. Such a Bayesian inference is now referred as multi-point calibration, in contrast with single-point calibrations where only data from a single scenario are used.

From a numerical point of view, the inference is done by drawing samples from the prior pdf and the likelihood. For this purpose, we use a Markov-Chain Monte-Carlo (MCMC) sampler, based on the Metropolis-Hastings algorithm, which is well suited to represent non-classical pdfs. Precisely we use the implementation made available through the widely used pymc 1 python library. Details of the Metropolis-Hastings MCMC algorithm can be found in the user guide.

Typically, a large number of samples is required to converge the posterior distributions. The results presented hereafter were obtained by running samples from 10 6 to 10 7 draws, with a thin factor between 10 and 20, and a burn-in of about 50%. To check that the convergence is reached, we use the z-score method proposed by Geweke [START_REF] Geweke | Evaluating the accuracy of sampling-based approaches to calculating posterior moments[END_REF] available in the pymc package. It consists of comparing the mean and the variance of segments from the beginning and the end of the Markov chain:

= + z E[ ] E[ ] Var[ ] Var[ ] beginning end
beginning end [START_REF] Nitschke | Model-form and predictive uncertainty quantification in linear aeroelasticity[END_REF] This z-score is theoretically distributed as standard normal variate. It is thus a statistical tool based on the hypothesis that

= E[ ] E[ ] beginning end .
The point is that if this z-score falls within 2 standard deviations of zero, we cannot reject this hypothesis with 95% chance to be right. Along with this Geweke z-score, we also check for convergence through the evolution of the means, the auto-correlations, the correlations, the traces and the pdfs.

= = = = f M M S S D D f D D M M S S f M M S S f D D M M S S f M M S S d f D D M M S S f M M S S f D D M M S S ( | , , ) ( | , , ) ( | , ) ( | , , ) ( | , ) ( | , , ) ( | , ) ( | , ) i i j j j i i j i i j j i i j i i j i j i i j i i j j i j ( 26
)
where the subscript i of θ means that the random vector of parameters depends on the model and M is a discrete random variable defined on a subset

= … M M { , , } I 1
of all possible models. Consider now a QoI Δ and a prediction scenario S

. By means of the law of total probabilities [START_REF] Draper | Assessment and propagation of model uncertainty[END_REF] we then write:

= = = = = = = = = = f S f S D D M M S S p M M D D S S p S S ( | , , , ) ( | , , , ) ( | , ) (
) i I j J j i j i j j j 1 1 (27)
where p represents the probability mass function (pmf) of a discrete random variable and:

• = = = f S D D M M S S ( | , ¯, , ) j i
j represents the distribution of Δ obtained by propagating in the code the posterior distribution of θ i calibrated under model i by using data from scenario j,

• we assume independence between D and S.

As a consequence, the mean and variance of Δ are:

= = = E S E S D M S p M D S p S [ | , , , ] [ | , , , ] ( | , ) ( ) i I j J j i j i j j j 1 1 (28)
where we omitted to mention M, D and S to simplify the notations. In the preceding equations:

• E S D S [ | , ¯, , ] j j is computed with = = E S D S E S D M S p M D S [ | , , , ] [ | , , , ] ( | , ) j j i I j i j i j j 1 ( 30
)
represents the mean of Δ by taking into account all the possible models M calibrated on the same scenario S j ,

• p M D S ( | ¯, )

( | , ) ( | ) ( | , ) ( , ) i j j j i j i j i I j i j i j 1 ( 31
)
where p(M i |S j ) is a prior user-defined pmf and p D M S ( ¯| , )

+ I i=1 J j=1 E[∆|S , D j , M i , S j ] -E[∆|S , D j , M, S j ] 2 p(M i |D j , S j )p(S j )    between-models, within-scenario variance + J j=1 E[∆|S , D j , M, S j ] -E[∆|S , Z, M, S] 2 p(S j )    between-scenarios variance ( 29
) = p D M S f D M S f M S d ( | , ) ( | , , ) ( | , ) j i j j i i j i i j (32)
The prior model mass function is chosen equal to: [START_REF] Rusak | Low-drag airfoils for transonic flow of dense gases[END_REF] meaning that all models have the same prior probability before observing data. Finally we chose the prior scenario mass function as in [START_REF] Edeling | Predictive rans simulations via bayesian modelscenario averaging[END_REF]:

= p M S I ( | ) 1 i j
= = = = = p S S E S D M S E S D S ( ) (a) [ | , , , ] [ | , , , ] (b)
j j i I j i j j j 1 /
1 / 1 2 j p j J j p 1 [START_REF] Cramer | Transonic flows of bethe-zel'dovich-thompson fluids[END_REF] where ϵ j in Eq. ( 34b) is a measure of the dispersion of different model predictions for scenario S′ using coefficients calibrated under scenario S j .

A special subclass of BMSA problems arises when there exists only one model. The BMSA equations are simplified and reduced to

and such problems will now be referred as Bayesian Scenario Averaging (BSA). In this case the scenario pmf (34a) cannot be used any more. We suggest instead the following formulation:

where the first term of the right hand side represents the error between the posterior mean output calibration quantity y and the reference data D ¯, j and the second term stands for the model inadequacy.

Calibration results

In this section we apply the statistical calibration framework described in Section 4.1 to infer on the input coefficients of two thermodynamic models. As anticipated in Section 3.2, the data used for the calibration are synthetic pressure measurements at 17 locations along the airfoil wall. Previous sensitivity analyses [START_REF] Merle | Bayesian quantification of thermodynamic uncertainties in dense gas flows[END_REF] show that, for the PRSV model, the most influential parameters are the acentric factor ω and the specific heat c v, ∞ , with n having an almost negligible effect. For MAH, only the critical temperature and pressure T c , p c have a significant influence on the model variance, while the effect of other parameters is rather negligible. Also note that the two parameters exhibit significant interactions (i.e. varying them jointly has a strong effect on the variance) for flow cases with shock waves, due to the strong nonlinear effects characterizing the shock region. Such interactions have a strong impact on the calibration results, leading to posterior correlations of the calibrated parameters. Specifically, in our previous work [START_REF] Merle | Bayesian quantification of thermodynamic uncertainties in dense gas flows[END_REF] a clear posterior correlation between ω and c v, ∞ was observed when calibrating the PRSV EOS. Similarly, a correlation between T c and p c was observed for MAH. For this reason, in both cases we chose to fix one of the correlated parameters. The acentric factor ω and the critical pressure p c are then removed from the uncertain parameter vector of the PRSV and MAH models, respectively, and kept fixed to their nominal values. For the remaining uncertain parameters and the hyperparameters we assumed large non-informative prior distributions, given in Table 3, which encompass the nominal values. The upper and lower bounds allow for a wide search of the parameter space, while avoiding physically inadmissible values (e.g., negative values or values leading to thermodynamic conditions lying in the two-phase regionsee Section 3.2).

Results for the PRSV model

The results of the calibrations of the PRSV EOS against different sets of data are reported in Table 4 in terms of average and standard deviation of the posterior distributions obtained for the various scenarios. The calibration against the three data sets simultaneously (multiplepoint calibration) is referred to as scenario "123" in the following. In all cases the physical parameter c v and the hyperparameter σ are well informed from the data.

Concerning the posterior distributions of the PRSV parameter c v, ∞ , two tendencies are observed. For calibration scenarios 1 and 123 the posterior mean of c v, ∞ is close to 90.0, i.e. about 20% higher than the nominal value. For scenarios 2 and 3 the mean is close to approximately 136. More generally, calibration scenarios 2 and 3 lead to similar posteriors for all the parameters and hyper parameters. This is mainly because the operating points 2 and 3 are closeby in the Clapeyron diagram and the corresponding pressure data are similar. In all cases, calibration leads to extremely peaked posteriors for c v, ∞ , whose support is much smaller than the prior distributions. This indicates that the computed posteriors are independent on the selected prior. The posterior mean of the hyper parameter σ is the smallest for scenario 2 (E [σ] ≈ 0.262) and the largest for point 1 (E[σ] ≈ 0.734). We recall that, in our statistical models for the likelihood based on Eq. ( 23) or Eq. ( 24)

Table 3

Prior distributions of physical and hyper parameters for the two EOS under investigation. σ is a measure of the ability of a model to capture the data: the more the model fails to fit the experimental points, the larger is σ. In other words, if the posterior distribution of the code output y exhibits a weak probability to observe the data, then σ increases to compensate for this discrepancy. In practice, due to model inaccuracy, the statistical model tends to capture the data by increasing the variability of the posterior predictive distribution of the QoI. This is done either through a large posterior variability of the physical parameter (in this case, c v ), or through the model-inadequacy term. This is clearly illustrated by Fig. 3, which reports the posterior predictive distributions (p.p.d.) for the pressure coefficient C p . These are compared to the nominal solution and to the data used for the calibration. The expectancy of the p.p.d. is found to be closer to the data for scenario 1 than for scenarios 2 and 3. However, because of the small standard deviation of c v, ∞ and hence of C p for this scenario, the reference is assigned a poor chance of occurrence and the model-form error term leads to a large variance of the p.p.d. In this figure, the error bars associated with the parametric uncertainties (not reported) are of the order of line thickness, and almost superposed with the average solution. For all cases, the p.p.d. are in better agreement with the reference than the deterministic prediction based on the nominal parameters.

In Table 5 we report the L 2 norm of residual error with respect to the data for the p.p.d. of C p . The norm is extended to the 17 chordwise locations. The table shows that 1) independently on the data used for the calibration, the average posterior predictions always improve over the nominal solution; 2) except for the p.p.d. for scenario 1, the multiple-point calibration provides higher residual errors than the p.p.d. based on individual calibrations, due to the fact that the multi-point calibration represents a tradeoff among the observed data for the various scenarios.

Figure 4 illustrates the impact of the calibrations on the predicted thermodynamic behaviour. Specifically, the figure shows posterior predictions for the saturation curve and the transition line in the Clapeyron diagram. Reference data obtained from the SW model are also reported for reference. It is worth noticing that, according to the SW the = 0 curve is located below the saturation curve, D5 being not a BZT fluid according to this reference model. However, the fundamental derivative is rather close to zero for a wide range of conditions in the saturated vapour region (see Fig. 4(e)). Concerning the calibrated model, several considerations are in order. First of all, the saturation curve is not affected by the calibration, since its location depends only on the thermal equation of state and not on caloric properties like c v . Secondly, the average posterior location of the transition line is above the nominal one in all cases, indicating that the calibrated model predicts, in average, a wider inversion zone than the baseline. As a consequence, for scenarios 1-3, the fundamental derivative tends to take throughout the flow lower values with respect to the nominal model. This leads in turn to a reduction in the maximum Mach number and, as a consequence, to predictions that match more closely the reference data than the nominal solution. Lastly, the posterior prediction of the transition line varies significantly according to the scenario used for the calibration, as a consequence of the inadequacy of PRSV in reproducing .

( ) and posterior E C C PRSV [ | , ]

p p j ref ,

(

). bars (

± E C C PRSV S C C PRSV [ | , ] [ | , ] p p j ref p p j ref , ,
): parametric uncertainty y ( ).

d p ref , 1
( ), p C C PRSV ( | , ) d p ref , 2
( ), p C C PRSV ( | , ) d p ref , 3
( ), p C C PRSV ( | , ) d p ref , 123
( ).
the fluid thermodynamic behavior correctly in the whole flow field. Interestingly, if the EOS is calibrated directly by using thermodynamic data (and more specifically pvT data) [START_REF] Merle | Bayesian estimates of thermodynamic uncertainties for flows governed by dense-gas equations of state[END_REF] the inversion zone tends to disappear (as in the SW model) but predictions of the aerodynamic quantities based on the updated EOS are found to be less accurate than the nominal model. Indeed, in this case the fundamental derivative takes overall higher values throughout the flow, leading to stronger shocks.

Finally, the ability of the calibrated model to predict unobserved flows is investigated by propagating the parameters calibrated against scenarios 1-3 and 123 through scenarios 4 and 5. The results are depicted in Figures 5 and6, along with validation data obtained by running the reference SW model and the nominal solutions. The L 2 -norms of the errors with respect to the pseudo-experiments are also reported in panel e of both figures. Predictions are based on the sole propagation of the posterior parameter distributions through the new scenario, characterized by a small variance. This leads to small uncertainty bounds of the output QoI, which do not always encompass the validation data. This makes us argue that the parametric uncertainty alone is not sufficient to capture the "truth", and that the model-form uncertainty should also be taken into account. Unfortunately, the model-inadequacy term η can be hardly extrapolated to a scenario different from the calibration ones. In almost all cases, the average predictions are closer to the reference solution than the nominal solution. However, the improvement is more or less important depending on the calibration scenario used to train the parameters. As expected from the posteriors of c v, ∞ , we observe similar distributions for scenarios 2 and 3 on one side, and for scenarios 1 and 123 on the other side, for both prediction points. The predicted mean distributions based on calibrations against scenarios 2 and 3 show a good adequacy with the reference, with posterior errors of about 0.370-0.377 and 0.143 compared to the nominal errors of 1.730 and 1.719 for the nominal solutions at points 4 and 5, respectively. In particular, the shock location is drastically improved and lies very close to the reference. On the other hand, the predicted shock locations based on calibration scenarios 1 and 123 remain very close to the nominal solutions. Note that the prediction scenarios 4 and 5 are rather close to calibration scenarios 2 and 3, which explains the good performance of the model when calibrated against these points. Finally, even if the posterior mean distributions of C p are in closer agreement with the reference than the nominal solution, the error bars do not allow to capture the reference data all along the airfoils. This is due to the small posterior variance of c v, ∞ and to the lack of a model-form uncertainty term. In other terms, accounting for the parametric uncertainty only leads to a severe underestimation of the solution variance, more notably in the shock region. This point will be discussed further in the following.

To complete the analysis, we finally consider the p.p.d. for a global performance parameter which was not directly informed from the data, i.e. the drag coefficient, given by

= C D V d 1 2 2 with = D p dS n i • airfoil
the pressure force component in the flow direction i. The results are reported in Fig. 7 for predictions at conditions 4 and 5, based on parameter posteriors calibrated on scenarios 1, 2, 3 and 123, respectively. Table 6 provides the means and the standard deviations for the various p.p.d.s, along with the value predicted by using the nominal values of the parameters and the pseudo-experimental reference.

Due to the linear relationship between C d and p, the drag coefficient follows similar trends as the pressure coefficient C p . However, due to the sensitivity of this integrated parameter to small variations of the shock location, C d distributions are more affected by the parametric uncertainty. For both prediction scenarios, the average C d is to within 3% and 4% of the reference value, in the best cases (prediction based on scenarios 2 and 3), which represents a considerable improvement over the nominal model results, characterized by relative errors of 21.18% and 26.88%, respectively. The prediction error increases to about 15% when using posteriors calibrated from scenarios 1 and 123. Despite the improvements observed for the prediction of mean values, it appears that the reference values are in general not captured by the posterior distributions of C d for any choice of the calibration scenario. In the best cases, the reference value of C d is close to the tail of the posterior distribution.

Results for the MAH model

In this section we discuss the calibration results for the MAH model. The posterior means and standard deviations of the stochastic parameters/hyperparameters T c , σ, α and β are given in Table 7 for the various calibration scenarios (single-and multi-point). Posterior distributions of the critical temperature T c are found to be rather similar

E C MAH [ | , ] p j ref , - - - 0.038 S C MAH [ | , ] p j ref , - - - 0.023
for all calibration scenarios. Indeed, the posterior means range from 565.746 up to 574.720. This are about 8% lower than the nominal value of 619.15. The posterior standard deviations are very small (less than approximately 0.3%) for all cases, showing that this parameter is very well informed from the data, but also that the calibration is very sensitive to the calibration scenario, since the different posteriors do no overlap. Concerning the hyperparameters, the posterior mean of σ ranges from 0.247 (for scenario 2) up to 0.837 (for scenario 123). These rather high values indicate that the model-form error plays a crucial role for capturing the data, thus compensating the small posterior variance of the physical parameter T c . The posterior distributions of the pressure coefficient for the calibration scenarios, reported in Fig. 8, show a good agreement with the data. For scenarios 1 and 2, the data are located to within one experimental standard deviation from the mean posterior prediction, included in the shock region. For point 3, the shock location is predicted less accurately and some points belong to weak probability regions of the posterior. Nevertheless, the mean posterior distribution represents a clear improvement with respect to the nominal solution. This is better seen by inspection of the L 2 norms of the residual error with respect to the data, given in Table 8. For all calibration scenarios and prediction points the error is reduced by approximately 20% (in the worst case) to approximately 90% (in the better case). Similarly to the PRSV model, the multiple data approach yields results close to the single data calibrations, with posterior values of T c close to those found for scenario 1. This seems to indicate that data from scenario 1 tend to dominate the calibration.

Posterior distributions of the fluid properties in the Clapeyron diagram are reported in Figures 9(a) to 9(d), along with the reference solution (Figure 9(e)). For the saturation curve, only the mean of the posterior distribution is reported, since the results are practically insensitive to T c . For the transition line, the variability is also very small and the error bars are drawn by applying a large multiplicative coefficient to the standard deviation. Contrary to the PRSV case, calibration of the MAH model leads to a significant shrinking of the = 0 curve with respect to the baseline. As a results, the posterior average area of the inversion zone is equal to zero for all calibration scenarios (Table 9(f)). Indeed, the nominal MAH model predicts a shock much weaker than the pseudo-experiment, and located upstream of the reference, i.e. an opposite situation with respect to the PRSV case. The calibration acts in the sense of increasing the shock strength by moving the thermodynamic parameter toward a range of values that lead to a smaller (or null) inversion zone and to overall higher values of Γ throughout the flow.

The posterior predictions of the C p are reported in Figures 10 and 11 for scenarios 4 and 5, respectively.

Due to the low sensitivity of the posterior distributions of T c to the calibration scenario, for both prediction scenarios similar posterior predictive distributions of the pressure coefficient are obtained for any choice of parameter distributions. The shock location is not perfectly captured but is significantly improved with respect to the baseline model. As already observed for the PRSV case, assimilating data from several scenarios at once does not really improve predictive accuracy over single-point calibrations. We also observe that the calibrated MAH model leads to less accurate predictions of the new scenarios 4 and 5 than PRSV, showing once more the low robustness to parametric uncertainties of this more complex model.

Finally, posterior distributions of the drag coefficient, reported in Fig. 12 and Table 9, exhibit a similar behaviour to the PRSV case. Calibration leads to a very significant improvement of C d predictions over the baseline model for both prediction scenarios, although the reference value is not captured by the posterior distributions (or is captured with very low probability), the best results being obtained when using posteriors calibrated on scenario 2.

Bayesian mixture model results

The results of Section 5 show that calibrating the parameters of thermodynamic models on one or more scenarios does not lead to robust predictions of new cases, even when these are not very far apart from the calibration scenarios. Not only the average prediction is not in close agreement with the available validation data but, more importantly, the predictive variance is strongly underestimated, since the sensitivity of a given model to calibration scenarios is not accounted for. Furthermore, the standard calibration approach does not provide any estimate of model-form errors for the predictive setting.

To overcome these drawbacks, in this section we carry our predictions of the unseen scenarios 4 and 5 by using Bayesian mixture models. As a first step, we incorporate into the models information coming from different calibration scenarios by averaging predictions based on different calibration data. This corresponds to the mixture model that we called BSA. The results are compared to those of the model calibrated simultaneously on scenarios 123. Afterwards, we average predictions of the PRSV and MAH models trained on various scenarios via the BMSA approach, and we assess the capability of the mixture model to capture validation data.

Accounting for scenario uncertainty: BSA results

The BSA mixture approach is first applied to the PRSV EOS. The posterior distributions obtained for calibration scenarios 1, 2 and 3 are propagated through a new scenario, and the posterior predictive distributions are averaged by using the a priori scenario p.m.f. (34a). A preliminary sensitivity study was carried out to investigate the influence of the exponent p. We recall that using = p 0 corresponds to as- signing the same probability to all calibration scenarios, whereas high values of p tend to assign a higher probability to the scenario for which the chosen error measure is the lowest. The results (not reported for brevity) show that the predictions tend to become independent of p for p ≳ 5. For such high values, the p.m.f. assigns a probability very close to one (90% or more) to the scenario with the lowest error measure for each case, so that the BSA prediction becomes equivalent to a singlepoint standard prediction. This is due to the fact that the error measure is based on the model performance for the calibration scenario, data for the prediction one being not available. Additionally, the predictive variance reduces to the one of single-point predictions, i.e. does no longer account for scenario uncertainty. Using low values of p (comprised between 0 and 2 according to the case) leads to more accurate predictions of the expected solution, while providing a more conservative estimate for scenario uncertainty. In the following we restrict our analyses to the case = p 1. The resulting prior scenario probabilities are reported in Table 10.

PRSV-based BSA predictions of the wall pressure coefficient for scenarios 4 and 5 are reported in Figure 13, along with the reference data C p ref ,{4,5} and the nominal PRSV solution C p Nom ,{4,5} . Predictions based on the multi-point calibration (scenario 123) are also reported for comparison. For both predictions, the average BSA solutions are close to the reference data, except in the vicinity of the shock. BSA not only outperforms the baseline model, but is also significantly more accurate than the multi-point model calibrated on the same datasets. Additionally, BSA now provides a sufficiently conservative estimate of the posterior predictive variance, which encompasses the pseudo-experimental data.

Results based on the MAH model are reported in Fig. 14. For this model all results are in less good agreement with the reference than for PRSV. In the shock region, the predicted solution error bars ( ± 1 standard deviation) do not encompass the reference data. Nevertheless, BSA provides once again more accurate results and a more realistic estimate of the predictive uncertainty than the standard calibration approach, as seen by inspection of the prediction errors (Table 11). This is also reflected by the predicted drag coefficient, reported in Table 15 for both PRSV and MAH models. The expected value of C d according to the BSA model is much closer to the reference, both with respect to the nominal model and the multi-point calibration. Specifically, the BSA prediction is to within approximately one standard deviation from the reference value, whereas with the standard approach the reference was not captured by the posterior predictive distribution.

Accounting for model-form uncertainty: BMSA results

The BSA approach allows to incorporate information about the sensitivity of model predictions to the calibration scenarios. The preceding results show the interest of mixing several scenarios to improve estimates of the posterior prediction of a given QoI, both in terms of expectancy and standard deviation. It was also seen that different models (PRSV and MAH) may lead to non overlapping posterior predictions. It is then necessary to incorporate in the prediction the notion of model-form uncertainty. This is done by considering the full BMSA formulation presented in Section 4.2. The latter is applied to predict scenarios 4 and 5 by mixing forecasts based on both PRSV and MAH models and the various calibration scenarios available. For the full BMSA formulation it is also possible to use the more general measure of the predictive error based on model agreement.

Table 13provides the posterior probabilities of the two models for each single-point calibration scenario, as well as the scenario probabilities. The Table also reports the posterior model probability for the multi-point case. In this case, only a single scenario is considered, which is assigned a probability equal to one, and BMSA reduces to BMA. Interestingly, while in the BMSA approach both model contribute to the posterior prediction, although the MAH model is found to provide a nonzero contribution only for scenario 2, only the PRSV model is picked by BMA applied to the multi-point calibrated models. Wall distributions of C p are reported in Fig. 15 for BMSA based on singlepoint scenarios and BMA based on the multi-point calibration. Using multiple scenarios and models not only provides a much more accurate prediction of the QoI, but also of the posterior variance, which is found to encompass the data for all cases. We observe that, contrary to the predictions based on the Gaussian model-inadequacy term (Figs. 3 and8), the predictive variance is now of the same order of the pseudoexperimental one in smooth flow regions and becomes larger in the shock region, where the predictions are very sensitive to the thermodynamic model and choice of coefficient. On the contrary, the BMA/ multi-point prediction underestimates both the parametric and the model-form uncertainty, leading to very inaccurate mean predictions in the shock region, as well as to a dramatic under-prediction of the posterior variance. The overall prediction errors for the two approaches and the two prediction scenarios are reported in Table 14.

The accurate prediction of the C p distribution along the wall and of the associated uncertainties enables a more reliable estimate of the drag coefficient (see

).

closer to the reference for both prediction scenarios. Additionally, the reference is to within approximately one standard deviation from the BMSA average prediction. On the contrary, not only the BMA/multi-point prediction is far apart from the reference, but it is also characterized by an unrealistically small posterior variance, i.e. the model assigns a very high probability to the wrong result.

Conclusion

A robust Bayesian prediction methodology for flows of dense gases subject to thermodynamic uncertainties has been described and assessed. As a first step, Bayesian inference was used to calibrate the input parameters of two popular equations of state (EOS), namely, the cubic Peng-Robinson-Stryjeck-Vera and the virial Martin-Hou equation of state. The inference was carried out by performing inverse propagation of a set of observed aerodynamic data (wall-pressure measurements) for transonic dense gas flows past an airfoil through a dense gas solver supplemented with one or the other EOS. Different data sets, corresponding to various calibration scenarios where used to calibrate the coefficients, allowing to investigate the sensitivity of the posterior parameter distributions to the observed data. Specifically, different scenarios were generated by changing the free-stream thermodynamic conditions of the incoming flow. The coefficients were calibrated both considering the calibration scenarios one by one or all of them simultaneously (single-point or multi-point calibrations). The latter approach was supposed to lead to more robust predictions, since the calibration is the result of a compromise on fitting the data available for different scenarios. Stochastic observational-error and model-inadequacy terms were introduced to construct the likelihood function. In all cases, the physical model parameters and the hyperparameters associated to the stochastic models were well informed by the data. The posterior parameter distributions were found to be rather sensitive to the calibration scenario, showing the lack of universality of the estimated parameters, which lose their physical meaning and become just tuning parameters, which are likely to incur overfitting problems. The parameter sensitivity to the calibration scenario is a result of model- form inadequacy. Despite calibration, these was found to be so large that neither model was able to fit the data perfectly well. However, in most cases these were captured by the tail of the model-inadequacy distribution. Then, the various posterior distributions were used to predict unobserved flow scenarios and quantities of interest (QoI) different from the ones used for calibration. Since the model-inadequacy term cannot be transposed to different flows or QoI, the predictions were carried out by just propagating the parameter posteriors through the dense-gas solver. The results were compared to the available validation data and to deterministic predictions based on the nominal parameter values. Although the predictions of the calibrated models are significantly more accurate than the baseline, they did not fully agree with the reference data. Additionally, the predictive posterior variance was found to be too small to encompass the data. Multi-point calibration did not really improve the predictions over single-point ones.

To better account for scenario uncertainty while improving the robustness of the predictions, several Bayesian mixture models were considered. When a single candidate model is available, a Bayesian Scenario Averaging (BSA) is carried out, i.e. a mixture of the predictions based posteriors calibrated on different scenarios, applied to a new case. Since the prediction scenario is independent on those used for calibration, the weights are prior scenario probabilities. These are constructed by assigning higher weights to scenarios for which the calibration fits well the observed data. Mixing scenarios provides a much more accurate prediction of the new case, both in terms of posterior expectancy and standard deviation. Specifically, BSA provides a much more realistic estimate of the posterior variance due to model sensitivity to the calibration scenario than just propagating a posterior distribution of the parameters inferred from all the scenarios at once. However, BSA based on different models still lead to non-overlapping predictions for the new cases, meaning that the prediction is still sensitive to the chosen model structure. A more general prediction model is then constructed by averaging not only on the scenarios but also on the ). For BMA the between-scenario variance is identically zero.

c

  c denotes critical-point quantities, v the specific volume, and: the critical compressibility factor and R the gas constant. The function α(T r ) is an adimensional relationship depending on the reduced temperature and the substance acentric factor ω:

  the solution of a cubic equation for Z c with only one relevant root ( = Z 0.3112 c

  T c (K) 619.15 p c (atm) 11.45 Z c 0.286 T e (K) 484.1 ω 0.6658 c ν, ∞ (T c )/R 76 n 0.5208 distribution of the pressure coefficient along the airfoil (Fig. 1(b)).

Fig. 1 .

 1 Fig. 1. Typical reference solution and location of the numerical pressure taps.

Fig. 2 .

 2 Fig. 2. Representation in the Clapeyron diagram of the operating points (scenarios for the free-stream thermodynamic conditions) used in the pseudo-experiments. SW reference ( ), PRSV ( ), MAH (). Red symbols connected by lines: calibration scenarios; green symbols: prediction scenarios. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

4. 2 . 1 .

 21 Bayesian model-Scenario averaging (BMSA)Let us now consider = … i I 1, , models, modelling the same physical problem, involved in = For model i under scenario j, we write the Bayesian rule as:

j

  is the posterior model probability, computed by applying again Bayes rule = = p M D S p D M S p M S p D M S p M S ( | , )

  ∆|S , D j , M i , S j ]p(M i |D j , S j )p(S j )

Fig. 4 .

 4 Fig. 4. Calibration results for PRSV. Posterior predictions for the saturation curve and the transition line = ( 0) . Panels (a) to (d): nominal saturation curve ( ), nominal transition line ( ), mean posterior saturation curve (+), mean posterior transition line ( × ) and posterior transition line error bars ( • ). Panel (e): reference data based on the SW model ( ), transition line ( ), = ( ) and = ( ). Panel (f): posterior average and standard deviation of the inversion zone area. For the nominal PRSV model, the area of inversion zone is equal to 6.429 10 4 .

Fig. 5 .Fig. 6 .

 56 Fig. 5. Prediction results for the calibrated PRSV model. Pressure coefficient for operating point 4: reference data C p ref ,4 ( ), nominal C p Nom ,4 . ( ) and posterior E C C PRSV [ | , ] p p j ref ,

Fig. 7 .

 7 Fig. 7. PRSV -Drag coefficient predictions. Reference C d ref ,{4,5} ( ), nominal C d Nom ,{4,5} . ( ), p C C PRSV ( | , )

Fig. 8 .

 8 Fig. 8. MAH model: posterior predictive distributions (p.p.d.) of the pressure coefficient based on various calibration scenarios. Reference C p j ref , ( ), nominal C p j Nom , .

Fig. 9 .

 9 Fig. 9. Calibration results for MAH. Posterior predictions for the saturation curve and the transition line = ( 0). Panels (a) to (d): nominal saturation curve ( ), nominal transition line ( ), mean posterior saturation curve (+), mean posterior transition line ( × ) and posterior transition line error bars ( • ). Panel (e): reference data based on the SW model ( ), transition line ( ), = 0.3 ( ) and = 0.5 ( ). For the nominal MAH model, the area of the inversion zone is equal to 1.391 10 2 .

Fig. 13 .(

 13 Fig. 13. PRSV -BSA prediction of the pressure coefficient. Reference C p ref ,{4,5} ( ), nominal C p Nom ,{4,5} . ( ) and posterior = E C M PRSV [ | ] p

Fig. 14 .

 14 Fig. 14. MAH -BSA prediction of the pressure coefficient. Reference C p ref ,{4,5} ( ), nominal C p Nom ,{4,5} . ( ) and posterior = E C M MAH [ | ] p

Fig. 15 .

 15 Fig. 15. BMSA predictions of the pressure coefficient compared to the results of BMA combined with multi-point calibration. Reference C p ref ,{4,5} ( ), BMSA/BMA posterior E[C p ]( ). Error bars: within-model within-scenario ( ), between-models within-scenario ( ), between-scenarios (for (a) and (c) only) ( ), total (). For BMA the between-scenario variance is identically zero.

Table 2

 2 Numerical values: specific volume v, density ρ end pressure p.

	Point	v/v c	ρ/ρ c	p/p c
	1	1.60	0.625	0.970
	2	2.00	0.500	0.950
	3	2.10	0.476	0.900
	4	1.90	0.526	0.940
	5	1.84	0.543	0.933

Table 4

 4 PRSV model: mean E and standard deviation S of the posterior distributions of the parameters for the various scenarios.

	Mean and Std. deviation	Scenario j
												1	2	3	123
	E c [	v	,				C | p j ref ,	,	PRSV	]	90.772	135.762	135.940	91.827
	S c [	v	,				C | p j , ref	,	PRSV	]	0.247	6.013	4.026	0.555
	E C [ | p j , ref	,	PRSV	]	0.734	0.262	0.447	0.379
	S C [ | p j , ref	,	PRSV	]	0.773	0.067	0.114	0.049
	E C [ |	ref	,	PRSV	]	50053.635	0.033	0.034	0.011
								p j ,
	S C [ |	ref		,	PRSV	]	28912.708	0.016	0.013	0.003
						p j ,
	E C [ |	ref		,	PRSV	]	-	-	-	0.123
							p j ,
	S C [ |	ref	,	PRSV	]	-	-	-	0.028
					p j ,

Table 5

 5 Calibration results for PRSV -L 2 -norm error of the residual error with respect to the reference C p data.

	Scenario	C p j ref ,	C p j Nom ,	.	2	C p j ref ,	E C C [ | p p j ref ,	,	PRSV	]	2	C p j ref ,	E C C [ | p p ref ,123 ,	PRSV	]	2
	j					(Single-point calibration)		(Multi-point calibration)
	1	1.273				0.130						0.115
	2	1.393				1.023						1.032
	3	1.698				1.274						1.720

Table 6

 6 PRSV -Predictions of the drag coefficient.

	Prediction	Calibration	Reference	Nominal	Posterior						Relative error (%)
	scenario k	scenario j	C d k ref ,	C d k Nom ,	.	E C C [ | d p j ref ,	,	PRSV	]	S C C PRSV [ | | , d p j ref	]	Nominal	Posterior
	4	1	1.223 10 1	1.482 10 1	1.414 10 1				1.075 10 4		21.18	15.62
	4	2	-	-		1.260 10 1				1.702 10 3		-	3.03
	4	3	-	-		1.260 10 1				1.055 10 3		-	3.03
	4	123	-	-		1.410 10 1				2.268 10 4		-	15.29
	5	1	1.079 10 1	1.369 10 1	1.265 10 1				1.617 10 4		26.88	17.24
	5	2	-	-		1.040 10 1				2.398 10 3		-	3.61
	5		-	-		1.039 10 1				1.471 10 3		-	3.71
	5	123	-	-		1.258 10 1				3.393 10 4		-	16.59

Table 7

 7 MAH model: mean E and standard deviation S of the posterior distributions of the parameters for the various scenarios.

	Mean and Std. deviation	Scenario j
										1	2	3	123
	E T C [ | c p j , ref	,	MAH	]	573.898	565.746	572.271	574.720
	S T C [ | c p j , ref	,		MAH	]	1.604	1.485	0.462	0.689
	E C [ | p j , ref		,		MAH	]	0.474	0.247	0.772	0.837
	S C [ | p j , ref	,	MAH	]	0.104	0.056	0.773	0.131
	E C [ |	ref	,	MAH	]	0.009	0.009	50182.526	0.017
		p j ,
	S C [ |	ref		,	MAH	]	0.004	0.005	28864.872	0.004
	p j ,

Table 8

 8 Calibration results for the MAH model -L 2 -norm error of the pressure coefficient with respect to reference C p data.

	Scenario j	C p j ref ,	C p j Nom ,	.	2	C p j ref ,	E C C [ | p p j ref ,	,	MAH	]	2	C p j ref ,	E C C [ | p p ref ,123 ,	MAH	]	2
						(Single data calibration)		(Multiple data calibration)
	1	2.255				0.429						0.316
	2	1.808				0.115						0.637
	3	2.289				1.761						1.785

Table 15 )

 15 . Indeed, BMSA predictions of C d are much Fig. 10. MAH. Prediction results for the calibrated MAH model. Pressure coefficient for operating point 4: reference C p

	ref ,4 ( ), nominal C p Nom ,4	. (	) and posterior

Table 9

 9 MAH -P redictions of the drag coefficient.

	Prediction	Calibration	Reference	Nominal	Posterior						Relative error (%)
	scenario k	scenario j	C d k ref ,	C d k Nom ,	.	E C C [ | d p j ref ,	,	PRSV	]	S C C PRSV [ | | d p j ref ,	]	Nominal	Posterior
	4	1	1.223 10 1	5.817 10 2	1.002 10 1				2.368 10 3		52.44	18.07
	4	2	-	-		1.126 10 1				2.541 10 3		-	7.93
	4	3	-	-		1.025 10 1				6.767 10 4		-	16.19
	4	123	-	-		9.892 10 2			1.027 10 3		-	19.12
	5	1	1.079 10 1	2.593 10 2	7.069 10 2			3.052 10 3		75.97	34.49
	5	2	-	-		8.679 10 2			3.492 10 3		-	19.56
	5	3	-	-		7.359 10 2			8.608 10 4		-	31.80
	5	123	-	-		6.908 10 2			1.285 10 3		-	35.98

Table 10

 10 BSA -Scenario probabilities for the PRSV and MAH model ( = p 1).

		PRSV	MAH
	p(S 1 )	0.460	0.260
	p(S 2 )	0.309	0.648
	p(S 3 )	0.231	0.093

Table 11 L

 11 2 -norm of the prediction error ( = p 1). Nominal= C C ,

								p ref ,{4,5}	p Nom ,{4,5} .	2 Multi-
	point= C p ref ,{4,5}	E C C [ | p p ref ,123 2 and BSA= C ] p ref ,{4,5}	E C M [ | ] p	2 .
	Scenario	PRSV			MAH
	j	Nominal	Multi-point	BSA	Nominal	Multi-point	BSA
	4	1.730	1.044	0.583	2.449		1.342	1.181
	5	1.719	1.697	0.787	2.649		1.905	1.704

Table 12

 12 BSA -posterior drag coefficient for = p 1.

	Scenario	Reference	Posterior PRSV		Posterior MAH
	j	C d j ref ,	E[C d |PRSV]	S[C d |PRSV]	E[C d |MAH]	S[C d |MAH]
	4	1.223 10 1	1.331 10 1	7.774 10 3	1.084 10 1	6.147 10 3
	5	1.079 10 1	1.143 10 1	1.078 10 2	8.139 10 2	8.038 10 3

Table 13

 13 BMSA -Posterior model probabilities and prior scenario probabilities.

	Scenario j	p(S k )		p(PRSV|S j )	p(MAH|S j )
		Scenario 4	Scenario 5
	1	0.291	0.266	1.000	0.000
	2	0.382	0.385	0.809	0.191
	3	0.327	0.349	1.000	0.000
	123	1.000	1.000	1.000	0.000

Table 14 L 2

 142

	-norm of the posterior prediction error: C p ref ,{4,5}	E C [ ] p	2 for BMA
	with multi-point calibration and BMSA.
	Scenario	BMA +multi-point		BMSA
	j	calibration
	4	1.044		0.382
	5	1.697		0.497

Table 15

 15 BMSA predictions of the drag coefficient compared to the results of BMA combined with multi-point calibration.

	Scenario	Reference	Posterior (BMSA)	Posterior (BMA+multi-point)
	j	C d j ref ,	E[C d ]	S[C d ]	E[C d ]	S[C d ]
	4	1.223 10 1	1.273 10 1	0.772 10 2	1.410 10 1	2.268 10 4
	5	1.079 10 1	1.059 10 1 1.069 10 2	1.258 10 1	3.393 10 4

competing models (Bayesian Model-Scenario Averaging, BMSA). In this case the models are weighted by their posterior probabilities (resulting from inference on a given scenario) while the scenarios are assigned a prior probability mass function based on the level of agreement of the models for the new case, rather than model fitness for the calibration case as in BSA. Although one of the competing models (namely, Martin-Hou) is always assigned low probability compared to the other one, still its contribution to the mixture allows not only to improve the average posterior prediction, but also provides a more accurate estimate of the variance. Comparisons with the results provided by the Bayesian model averaging (BMA) of the two EOS with coefficients calibrated on several scenarios simultaneously clearly show the importance of scenario averaging for achieving a correct prediction of the variance. Indeed, scenario uncertainty is much more significant than model-form uncertainty in the present application, due to the dominant role taken over by one model over the other. In conclusion, BMSA appear to be an effective tool for the prediction of dense gas flows with quantified thermodynamic uncertainty.