Uncertainty Quantification of GEKO Model Coefficients on Compressible Flows

In the present work, supersonic flows over an axisymmetric base and a 24-deg compression ramp are investigated using the generalized k-ω (GEKO) model implemented in the commercial software, ANSYS FLUENT. GEKO is a two-equation model based on the k-ω formulation, and some specified model coefficients can be tuned depending on the flow characteristics. Uncertainty quantification (UQ) analysis is incorporated to quantify the uncertainty of the model coefficients and to calibrate the coefficients. The Latin hypercube sampling (LHS) method is used for sampling independent input parameters as a uniform distribution. A metamodel is constructed based on general polynomial chaos expansion (gPCE) using ordinary least squares (OLS). The influential coefficient closure is obtained by using Sobol indices. The affine invariant ensemble algorithm (AIES) is selected to characterize the posterior distribution via Markov chain Monte Carlo sampling. Calibrated model coefficients are extracted from posterior distributions obtained through Bayesian inference, which is based on the point-collocation nonintrusive polynomial chaos (NIPC) method. Results obtained through calibrated model coefficients by Bayesian inference show superior prediction with available experimental measurements than those from original model ones.


Introduction
Recently, uncertainty quantification (UQ) has been applied to computational fluid dynamics (CFD) by improvement of computing performance and reduction of computational cost. UQ is used to estimate not only the mean value and standard deviation but also the probability distributions of QoIs (Quantities of Interest) by considering the probability distributions of input variables such as boundary conditions or applied model coefficients. The Monte Carlo (MC) technique has been the general method used to compute output variables by sampling the random inputs with a specific distribution. Its basic concept is relatively simple but too expensive to achieve reasonable convergence, which means that application to CFD is difficult. The generalized polynomial chaos (gPC) method that was recently proposed by Xiu [1] requires much fewer samples than the MC technique so that studies using UQ can be performed efficiently.
Uncertainties in CFD can be divided into aleatory uncertainty and epistemic uncertainty. Aleatory uncertainty relates to the random components contained in a real system; from the perspective of CFD, examples would include uncertainties in geometry and boundary conditions such as velocity and pressure. On the other hand, epistemic uncertainty occurs due to a lack of knowledge regarding physical phenomena and modeling hypotheses. Thus, uncertainties caused by the process of model formulation and estimation of model coefficients, such as those used in modeling turbulence or multiphase flow, are typical examples of epistemic uncertainty.
The methodology of the UQ can be classified into two approaches, the forward problem and the inverse problem, depending on how the random inputs are treated [2]. In the case of a forward problem, the distribution of random inputs is given with a fixed mean and interval, and then one needs to compute the stochastic outputs with mean and deviation values based on the selected UQ methodology. However, in the inverse problem, the Bayesian inference is adopted to estimate the distribution of random inputs using observation data such as experimental results or synthetic data. Through the process of model calibration, the distribution of the random inputs is updated, and then QoIs are computed statistically; this process is called model prediction.
Schaefer et al. [3] applied point-collocation nonintrusive polynomial chaos (NIPC) to various turbulence models to study the effect of the random inputs, including aleatory and epistemic uncertainties, on the aerodynamic performance coefficients. They assumed model coefficients of the Spalart-Allmaras (SA) model and k-ω family models as a uniform distribution and obtained the statistical distributions of pressure and aerodynamic coefficients over an RAE 2822 airfoil and an axisymmetric bump. Additionally, influential coefficients among the model coefficients of the corresponding models are investigated through Sobol indices.
For high-speed flows, Huan et al. [4] conducted a global sensitivity analysis to identify effective input parameters, which are inflow, fuel inflow, and wall boundary conditions, and turbulence model parameters in the CFD of a scramjet system by using UQ methods. Additionally, the estimated model errors for models of different fidelity were investigated. Burt and Josyula [5] considered aleatory and epistemic uncertainties in sensitivity analysis/UQ calculations. Global sensitivity analysis and UQ are integrated with a direct simulation, which is a Monte Carlo gas flow simulation code for a hypersonic double-cone flow.
Recently, the generalized k-ω (GEKO) model, which is a two-equation model based on the k-ω model formulation, was proposed by Menter et al. [6]. The GEKO model itself includes model coefficients that can be controlled by users depending on the flow type, which means that the coefficients are able to be tuned in a variety of flows. Figat and Piatkowska applied the GEKO model to investigate the interaction between the propeller and the fuselage of an aircraft [7]. Brezgin and Aronson adjusted the value of C SEP for the sensitivity regarding adverse pressure gradients and spreading rates, which were overpredicted by conventional models (standard k-ω, SST k-ω) [8]. Additionally, Li et al. [9] compared the GEKO model with other turbulence models to investigate the influence of wall conduction on the heat transfer of supercritical n-decane in active regenerative cooling channels.
In the present work, supersonic flows over an axisymmetric base and 24-deg compression ramp are investigated. In the simulation of a 24-deg compression ramp, predicting shock wave/boundary layer interaction (SWBLI) accurately is an important factor in the design of high-speed flight vehicles. On the exterior flow over the aircraft, SWBLI can cause loss of control, some peaks on surface thermal loading. In internal flows, it can enhance distortion and pressure losses and can even cause catastrophic events leading to engine unstart. Rizzetta [10] investigated supersonic flow over a 24-deg compression ramp by using explicit algebraic Reynolds stress models and k-ε models, but the model that gave the best overall results was the Speziale-Abid [11] realizable k-ε model. Gerolymos et al. [12] also used Reynolds stress models and concluded that the WNF-LSS RSM model gave the best prediction of skin friction at reattachment and in the relaxation region. In terms of experimental studies, Settles et al. [13] tested a 24-deg compression ramp and men-tioned that two-dimensionality can be considered because the three-dimensional perturbations of the flow are minor. However, aerodynamic fences are required to achieve twodimensionality. Horstman et al., Dolling and Murphy, and Settles and Dodson also experimentally studied compression ramps [14][15][16].
In the case of base flow, the accurate prediction of drag has been a challenge problem in CFD. A supersonic body experiences major drag from skin friction drag, wave drag, and pressure drag. If the drag is able to be controlled, stability and control of vehicles can be enhanced. As for an experiment with supersonic base flow, Herrin and Dutton [17] experimented on a supersonic axisymmetric base flow at M ∞ = 2:46. Many researchers have tried to analyze the supersonic base flow by using CFD. Forsythe et al. [18] used Reynolds-averaged Navier-Stokes (RANS) and detached eddy simulation (DES) with a compressibility correction to analyze supersonic base flow. Both models using the compressibility correction showed improved pressure distribution compared to the models that did not use the correction. In general, it is known that large eddy simulation and hybrid RANS/LES simulation are able to predict the pressure distribution and recirculation zone better than traditional RANS, but when the computational time and cost are considered, RANS is still the most widely used method in engineering.
In the present work, UQ with Bayesian inference is applied to the GEKO turbulence model to estimate the optimized model coefficients and QoIs in supersonic flows over an axisymmetric base and 24-deg compression ramp. The present work demonstrates Bayesian uncertainty quantification of GEKO turbulence model coefficients for the first time, to the best of our knowledge. The point-collocation NIPC technique [1,19], which is one of the typical nonintrusive methods in UQ, is adopted for the UQ framework. Three model coefficients are considered the random input variables. As for QoIs, the pressure coefficient along the base    [21] to solve compressible Navier-Stokes equations and the energy equation. The advection upstream splitting method (AUSM+) scheme is selected because it provides exact resolution of contact and shock discontinuities, preserves positivity of scalar quantities, and gives less dissipation than the Roe scheme does. As for spatial discretization, the thirdorder monotone upstream-centered (MUSCL) scheme is incorporated to improve spatial accuracy by reducing numerical diffusion. The Courant-Friedrichs-Lewy (CFL) number is set to gradually increase from 0.25 to 5 depending on the convergence circumstance. To compensate for the defect predicting low pressure on the base, a compressibility correction lowering the turbulent eddy viscosity production is applied to the computational models [18,[22][23][24]. On the other hand, for the 24-deg compression ramp, a compressibility correction is not used due to degraded prediction in near-wall regions [10,12,[25][26][27][28]. Results obtained with the compressibility correction are denoted by "-CC" in the manuscript. As men-tioned in the previous section, the GEKO model is adopted in the present work. Other models including baseline (BSL) k-ω model, shear stress transport (SST) model, and k-ε realizable model in FLUENT [21] are also used for comparison.

Turbulence
Model. The characteristic of the GEKO model is that the model coefficients can be controlled to tune the model to various flow scenarios. In this present work, three model coefficients are considered for compressible flows.
The GEKO model formulation is as follows: The coefficients of the GEKO model are implemented through the functions ðF 1 , F 2 , F 3 Þ, which can be controlled to achieve different goals in a variety of parts of the computational domain. Although the details of formulations ðF 1 , F 2 , F 3 Þ, including values of C SEP , C NW , and C JET , are not released, the range of the coefficient values and characteristics of each coefficient are given below: (1) C SEP . Main parameter for adjusting separation prediction for boundary layers. Affects all flows. Increasing C SEP reduces eddy viscosity, leading to greater sensitivity to adverse pressure gradients for boundary layers and to lower spreading rates for free shear flows.
(2) C NW . Affects mostly the inner part of wall boundary layers (no impact on free shear flows). Increasing C NW leads to higher wall shear stress and wall heat transfer rates in nonequilibrium flows.
(3) C MIX . Affects only free shear flows (boundary layer shielded due to function F blend , which is discussed later). Increasing C MIX increases spreading rates of free shear flows. For each value of C SEP , an optimal value of C MIX exists, which maintains optimal free shear flow. This value is given by the correlation C MIX = C MixCor , which is the default.
(4) C JET . Is active in a submodel of C MIX (no impact when C MIX is equal to 0). Affects mostly jet flows. Increasing C JET while C MIX is active decreases the spreading rate for jets.
Allows adjustments to the spreading rate of jet flows while maintaining the spreading rate of the mixing layer.
(5) C CORNER . Nonlinear stress-strain term to account for secondary flows in corners (e.g., wing-body junctions).
(6) C CURV . An existing model for curvature correction, which can be combined with the GEKO model.
All coefficients can be accessed globally or locally by using user-defined functions (UDFs), which allows a global or zonal model optimization.
The coefficients C MIX and C JET are designed for free shear flows, whereas C SEP and C NW have an effect on boundary layers. To avoid any influence of C MIX and C JET on boundary layers, a blending function is incorporated, which deactivates C MIX and C JET in the boundary layer. The blending function is given by This function activates following the shear flow parameters: There are two important aspects to F Blend . Firstly, inside boundary layers, the function F Blend is equal to 1, which means that the function F Free becomes equal to 0. Secondly,   International Journal of Aerospace Engineering the parameter C JET is a subparameter of C MIX . As mentioned above, it only affects the simulation in cases where C MIX is not equal to 0. The Min and Max values of C MIX are suggestions from [8]. There might be situations where values lower than Min (0.5) or higher than Max (1.0) would be appropriate for specific flows. To avoid negative effects on free mixing layers due to changes in C SEP , however, the use of C MixCor (Equation (3)) is recommended. C CORNER has an effect on the secondary flows at the corner in the 3D domain, and C CURV affects the cylinder geometry as a cyclone for the curvature correction. Hence, in the present work, model coefficients (C SEP , C NW , and C JET ) are varied within the ranges given in Table 1.

Point-Collocation Nonintrusive Polynomial Chaos (NIPC).
In this current work, point-collocation NIPC, which was proposed by Hosder et al. [29], is applied. The pointcollocation NIPC technique requires the minimum number of random input variables calculated by Equation (6) consisting of the polynomial order (p), the number of random input variables (n), and the oversampling ratio (n p ). The minimum amount of random data is extracted by using the Latin hypercube sampling (LHS) method. Extracted random input variables are used for Equation (7).
Here, α n ðXÞ are the polynomial chaos coefficients and Ψ n is an element of an orthogonal family, and Y represents the output corresponding to a given random input. Note that Y and Ψ n in the above equation correspond to given values, and α n is obtained by applying the Gram-Schmidt technique to the following matrix: The gPC technique utilizes the properties of a base function orthogonal. The basis functions of the gPC used according to the distribution of random input variables are shown in Table 2. In the present work, the Legendre basis function is used due to the uniform distribution.
The least square minimization method is adopted, which has the advantage that an arbitrary number of points can be used to calculate the coefficients, as long as they are a representative sample of the random input vector. Least square minimization is a method that minimizes the truncation error. Truncated PCE and a residual are expressed as follows: where P − 1 is the order of PCE, ε P is the truncation error, y α = fy 0 , ⋯, y P−1 g T is a vector containing the coefficients, and ΨðxÞ = fΨ 0 ðxÞ, ⋯, Ψ P−1 ðxÞg T is the matrix that assembles the values of all the orthonormal polynomials in X.
The least square minimization can be expressed as follows: Equation (10) is solved by ordinary least squares (OLS). The ordinary least square solution iŝ      International Journal of Aerospace Engineering Y is the model response corresponding to a random input vector.

Leave-One-Out Cross-Validation
Error ðϵ LOO Þ. The method of the leave-one-out cross-validation error ðϵ LOO Þ is adopted to evaluate the error of the constructed PCE. The ϵ LOO is designed to overcome overfitting limitations by using cross-validation. When the least square minimization is used for calculating the coefficients, the formulation to calculate ϵ LOO is where h i is the ith component of the vector, which is given below:  EXP [16] EXP [30] (c) EXP [16] EXP [30] (d)  In Equation (14), A is the experimental matrix. Mðx ðiÞ Þ is the model response and M PC ðx ðiÞ Þ is the PCE result. Here, b μ Y is the mean of the model response.

Bayesian Inference.
Bayesian inference is a process of updating the probability of θ based on the frequency at which datum x is observed after assuming a prior distribution of the random variable θ. Unlike the frequentist methodology, which estimates conditional probabilities based on a large number of samples, Bayesian inference estimates the posterior probability of a random variable with a small number of samples. The basic formula for Bayesian inference is shown below. Here, πðθÞ is the prior distribution of the input parameters. This could be set up heuristically or could be set up through a predetermined tendency. Similarly, πðx | θÞ is the likelihood. The likelihood measures the suitable statistical model with respect to given data that consider observations and prior distribution. The factor πðxÞ is used to normalize the right-hand side. As a result, the posterior distribution π ðθ | xÞ is finally calculated.
In the case of the independent observations given, the likelihood function is modeled as follows: If the independent observations are used, the likelihood function is expressed as a form of multiplication (Equation (16)). In Equation (18), Y i is an observation, M is a computational forward model, and x is a set of input parameters. The discrepancy term ∑ represents the effects of the measurement error, which is obtained from experimental references [4,[25][26][27][28][29].

Computational Setup
4.1. Geometry. The geometry of the base was selected from Forsythe et al. [18]. The cylinder length was set to 8R, where R is the base radius (31.75 mm). This cylinder length is determined to match the experimental momentum thickness [18]. The outlet is located 10R downstream from the base, while the far-field boundary is at 4:15R from the axis of symmetry.
In the case of the 24-deg compression ramp, the geometry was adopted from Gerolymos et al. [12]. A no-slip wall is placed on the top and bottom. The upper wall is located 0.2 m away from the bottom wall of the downstream section. The inlet is 0.2 m from the beginning of the 24-deg ramp, which is 0.15 m long. Various frames were used in the compression ramp experimental setups [13,30] for the profile measurements. The details of frames and specific locations where experimental data were measured are thoroughly explained in [12,13,30].

Grid.
A structured grid is generated by using Pointwise [31]. To resolve the turbulent boundary layer near the wall and satisfy the wall unit ðy + Þ within 1.0, which is averaged in the computational domain, the first cell height is set to 1 e − 06 (m). For the grid sensitivity test, a total of five types of grid were generated, which is discussed in detail in the result section. The total number of cells was 37,370 (100 × 75 upstream and 149 × 205 downstream from the base), while 180,000 cells (600 × 300) were generated for the compression ramp case. The grids are shown in Figures 1  and 2. . Adiabatic wall conditions were adopted at the wall.
In the case of the compression ramp, profiles obtained from the 2D channel simulation were specified for all dependent variables (e.g., k, ε, and ω) at the upstream inlet boundary of the computational domain. Supersonic inflow condition M ∞ = 2:85 was computed at the inlet boundary (a unit Reynolds number of Re = 63 × 10 6 m −1 ). The wall temperature was fixed to T w = 258:8 K, in accordance with measurements [13,14]. Extrapolation was applied at the downstream outlet boundary, and test conditions are shown in Table 3.

Grid Test.
The grid sensitivity was investigated using the grid convergence index (GCI) using five types of grid (tiny, coarse, medium, fine, and extra fine), which are characterized in Tables 4 and 5. There are several GCI methods, but simplified least square GCI (SLS-GCI), which was developed by Eça and Hoekstra [32], is used. The results obtained by SLS-GCI can be explained using three components with respect to QoIs:   (2) U e . The local error of estimation meaning the difference between the numerical and estimated results. To evaluate the grid test, RP and C Pbase (Equation (19)) are used for the axisymmetric base case. On the 24-deg compression ramp, RP and SP are considered. These are QoIs, which are used not only in the grid test but also for UQ. Table 6 shows U g and U e values. The values in the brackets show U e , and both are very close to zero, which means the computation is within the asymptotic range and the grids are reasonably generated. Table 7 and Figure 3 show extrapolation values with respect to the QoIs. The values inside the brackets in Table 7 mean the error (%) between each extrapolation value and the result of the medium mesh. Combination ð1, 3, 5Þ has the biggest spacing ratio, but the difference from the result of the medium mesh is less than 5% (Max 4.5%). Hence, the medium mesh is selected for the computational domain. Figure 4 shows the results of all grids.

Axisymmetric Base.
In the case of the base flow, the velocity profile was obtained 1 mm upstream of the base corner. Figure 5 shows the velocity profiles, which are compared with the profile obtained from the experimental data. Whether compressibility correction is used or not, the profiles of all models appear in a similar shape. Overall, all models predicted the boundary layer profile reasonably well under the current mesh. Pressure distributions along the base surface are shown in Figure 6. The base pressure distributions obtained without the compressibility correction are quite a bit lower than the experimental data [17]. On the other hand, the models including the compressibility correction show higher base pressure distributions but have greater variations along the base surface. Qualitative flat pressure distribution along the base surface is not predicted at all by RANS models [18, 22-24, 33, 34]. The reverse flow near the wake axis stagnates at the center of the base surface where the relatively high pressure appears (see Figures 6 and 7).
One of the QoIs for axisymmetric base flow is the shear layer RP. Figure 7 shows the axial velocity along the centerline of the axis. It is clearly shown that results from the models without compressibility correction match the experimental data well. If a smaller recirculation region is predicted, as in the models without compressibility correction, it causes the flow to turn sharply behind the base, leading to a more enhanced expansion wave and a reduction in pressure. A small separation region, therefore, causes larger pressure drag than a large separated region. Models with compressibility correction overpredicted not only the RP but also the peak reverse velocity. When the reverse flow is overpredicted, the flow is more accelerated outward along the base radius, which can cause a pressure reduction. The failure of RANS models to predict a flat pressure distribution along the base surface due to the large turbulent eddy viscosity is a well-known issue [34].

10
International Journal of Aerospace Engineering The contours of the Mach number and turbulent viscosity ratio are shown in Figure 8. The upper and lower parts of the figure correspond to GEKO-CC and GEKO, respectively. It is clear that the GEKO-CC model predicts RP further from the base than the GEKO model. As mentioned in the introduction, compressibility correction reduces turbulent eddy viscosity, which can lead to higher pressure levels, but it also increases the reattachment length. Thus, the size and length of the recirculation region behind the base are influenced by overpredicting or underpredicting the turbulent eddy viscosity. This factor can also be controlled by the C SEP coefficient of GEKO model, and it is indirectly indicated through turbulent viscosity ratio contours (Figure 8) how big an influence C SEP could have. The relationship between overprediction and/or underprediction of the RP and compressibility correction has also been investigated by many researchers [18, 22-24, 33, 34].

24-deg Compression Ramp.
For the 24-deg compression ramp case, the profile including all dependent variables is computed at the inlet. The boundary layer thickness (δ), displacement thickness (δ * ), and momentum thickness (θ) of the results obtained from the models (BSL, SST, and realizable models) are shown in Table 8. The profiles are extracted at a position where the profiles were measured in the experiment [30]. The SST model predicts a relatively lower value of δ than the other models do. Figure 9 shows nondimensional pressure distributions along the bottom surface of the computational domain. To compare computational results, experimental data obtained from [13,15,16] were used. (The experimental data was found in [16], but the exact source was not given [12].) The result of the GEKO model is quite similar to the result obtained from the SST model, as indicated [6]. Both models overpredicted the length of the reversed flow region, which means that a separation shock occurs too early. The models underpredicted even lower pressure compared to the experimental data after the corner. The BSL and realizable models performed the best in the pressure distribution. The initial pressure rise of the results obtained from both models matches the measured rise. This means that the separation shock predicted from numerical flow fields occurred in the same location as observed in the experimental flow fields.   Figure 10 shows the skin friction on the bottom surface of the compression ramp domain. As mentioned before, it is seen that the GEKO and SST models both overpredicted the separation length. Further, the GEKO and SST models predicted quite low skin friction in the relaxation region and more delayed RP than measured. In contrast with the GEKO and SST models, again, the BSL and realizable models show better agreement when it comes to the separation

12
International Journal of Aerospace Engineering shock, SP, and skin friction distribution after the RP. The BSL and realizable models accurately predicted the point that skin friction started decreasing, which means the separation shock location is predicted quite well. Looking at it more closely, the BSL model predicted the separation shock location slightly more accurately than the realizable model did, but the realizable model predicted the SP a little bit closer to the experimental data. As for the RP, all models failed to predict the RP perfectly. Having said that, however, the BSL and realizable models predicted the RP much more closely than the GEKO and SST models did, which means that the BSL and realizable models predicted much shorter separation lengths. The Mach number contours and static pressure are shown in Figures 11 and 12, respectively. In the Mach number contours, all models predict separation shock and reattachment shock, but the GEKO and SST models predict much larger separation bubbles than the BSL and realizable models do, as seen above in Figures 9 and 10. Additionally, the GEKO and SST models underpredict pressure compared to the realizable and BSL models in the relaxation region.

Uncertainty Quantification
5.3.1. ϵ LOO Sensitivity Analysis. The random input variables in the present work are C SEP , C NW , and C JET . The LHS method was used for sampling these input variables, assuming a uniform distribution. The selection of the interval of random inputs is critical in the UQ process and has many effects on the distribution of QoIs and the posterior one of RVs. In the present work, the interval of each model coefficient was set with the interval of coefficients of the GEKO model proposed by Menter et al. [6], which is shown in  13 International Journal of Aerospace Engineering of polynomial chaos and oversampling rate are set (Equation (6)). Once the total samples are computed in a deterministic solver and calculated, a metamodel (surrogate model) can be constructed by using the results obtained from the deterministic solver. One of the several standards validating the robustness of the constructed metamodel is ϵ LOO .
Even though the same number of samples was used, ϵ LOO of the axisymmetric base case is much higher than that of the 24-deg compression ramp case, which means that the axisymmetric base case is less robust than the 24-deg compression ramp case. To investigate what components affect ϵ LOO , we carried out a ϵ LOO sensitivity analysis. As mentioned above, the number of samples depends on the order of polynomial chaos (p) and oversampling rate ðn p Þ. The left side of Table 9 presents the number of samples required for each combination of p and n p . The right side of Table 9 shows the ϵ LOO corresponding to each QoI. ϵ LOO decreases only when the n p increases, which means increasing the n p could be a better choice if the data is limited or computational cost is expensive. However, an appropriate order for the polynomial chaos (p) is also crucial when it comes to Bayesian inference. Schaefer et al. [3] set n p = 2 and mentioned that most of the turbulence uncertainty can be quantified by using n p = 2. However, in the present work considering Bayesian inference, the oversampling rate n p = 2 was not sufficient to construct the surrogate model. In this study, n p is set to 3, and p is selected considering QoI distribution. Figure 13 shows the distributions of C P base and RP corresponding to the polynomial order, p. All distributions look quite similar to each other. In order to decide the specific order for the surrogate model and Bayesian inference, the 5th order is selected for the base flow. In the case of ramp flow, the 3rd order is chosen due to the limited data.

Global Sensitivity Analysis.
A total of 168 samples for base flow and 60 samples for ramp flow are used to construct the metamodel and to analyze which model closure coefficient is dominant in each flow. Table 10 shows Sobol indices of closure coefficients for QoIs of each test case. The largest contributors to uncertainty in each flow case are typed in italic (closure coefficients with Sobol indices of less than 3:0 × 10 −2 were not considered significant for reduced dimensionality analyses [3]). It is apparent that C SEP is the most dominant closure coefficient. As Menter et al. [6] noted, C SEP is the main parameter for adjusting separation prediction and has effects on the eddy viscosity for boundary layers and the spreading rate for free shear flows. The next dominant closure coefficient is C NW , which affects wall shear stress. Sobol indices of C NW in the 24-deg compression ramp case are higher than those in the axisymmetric base case. This could be presumed to be affected by the fact that the SWBLI phenomena occurred on the surface of the compression ramp. C JET has not made as much contribution to either case as the other closure coefficients have.

Bayesian Inference.
Based on the metamodel constructed using OLS, a likelihood function was calculated using the observation data. In the present work, experimental data of QoIs such as the pressure coefficient, RP in the base flow, and SP and RP in the ramp flow were applied for the observation data. Experimental data that is assumed to have a Gaussian distribution including an error of the experimental measurement is considered for the observation data.   The common finding from these four figures is that the prior distribution, which was assumed as a uniform distribution, produces a posterior distribution that does not show a uniform distribution any longer. This means that the results composed with the input variables assumed as a uniform distribution are not uniform distributions, and also, each input variable presents a specific distribution by the observation data, which were used to calculate the likelihood function. However, each model coefficient shows a different posterior distribution and different correlated values (red lines) depending on the selection of QoI, even though in the same flow, axisymmetric base flow, or compression ramp one. The correlated coefficients were extracted from the posterior distributions in Table 11, which shows the correlated coefficients according to QoIs in both the base and ramp flows. The correlated coefficients of C SEP are predicted as 2.290 for C P base and 0.723 for RP in the same axisymmetric base flow. When the minimum and maximum values of 0.7 and 2.5 are considered, two correlated values are located near the minimum or maximum values for different QoIs, C P base and RP. This means that if a better prediction of C P base is required, C SEP should be set to 2.290, and if the RP is predicted well, C SEP should be set to 0.723. Further, C NW , the second dominant parameter, is correlated differently and with a different sign. However, in the compression ramp flow, the values of C SEP are predicted similarly: 0.864 for SP and 0.713 for RP. When the default value of C SEP (1.75) is  Figure 18 present the results with the coefficients calibrated in terms of C P base and RP, respectively. Similarly, UQ-SP and UQ-RP in Figure 19 adopted the correlated coefficients for SP and RP, respectively. As shown in Figure 18, the pressure distribution of UQ-C P base shows the best agreement with the experimental data, with improvement over the results of the original GEKO-CC model and the UR-RP model since three model coefficients, C SEP , C NW , and C JET , are calibrated using C P base observation data for Bayesian inference. However, the RP is predicted best by UQ-RP, whose coefficients are calibrated with the RP, even though the lowest peak of the reverse velocity is overpredicted in UQ-RP relative to the others. As for the 24-deg compression ramp case in Figure 19, UQ-SP, which is calibrated with the SP data, shows better agreement in the pressure distribution and skin friction than the original GEKO model. UQ-RP predicts a delayed SP but more accurate RP relative to the other results. Tables 12 and 13 show the values of QoIs and the error between computational results and experimental data for each case. The values in the brackets indicate the percent error, and UQ-C P base shows a 4% error for C P base , which is an improvement of 2.5% over GEKO-CC. UQ-RP has a 7.5% error with respect to RP, whereas GEKO-CC shows a 25% error. In the compression ramp case, UQ-SP shows better agreement than GEKO-CC and UQ-RP, with a 3% error in the prediction of the SP. However, the error in RP is still high. UQ-RP is able to predict both SP and RP within 40%.

Conclusion
In the present study, the GEKO model was investigated to quantify the uncertainty of its coefficients for supersonic flows over an axisymmetric base and a 24-deg compression ramp. For the axisymmetric base flow, compressibility cor-rection aided the models to predict the pressure along the base surface more accurately than the models without any correction. However, using the compressibility correction increased the radial variation of the pressure due to the increased centerline velocity, which also increased the distance to the shear layer RP.
As for the 24-deg compression ramp case, in which the compressibility correction was not applied, the GEKO model failed to predict the separation shock point and separation length. It even underpredicted pressure and skin friction in the relaxation region.
To overcome these issues and quantify the uncertainties, PC-NIPC and Bayesian inference were applied to the coefficients of the GEKO model, which can be tuned for specific flows. It was found that the proper value of the oversampling rate in Bayesian inference with PC-NIPC might be 3 than 2 in the forward problem with PC-NIPC. The LHS method was used for sampling random input variables with a uniform distribution. Based on the surrogate model constructed by OLS, the dominant coefficient was able to be figured through Sobol indices. C SEP was the most influential coefficient in both test cases. Observation data made of Gaussian distribution considering experimental measurement errors was used to calculate the likelihood function. The MCMC was calculated by using AIES for better convergence in cases where a strong correlation between parameters is expected. It was found that the calibrated coefficients were computed differently depending on the flow type and QoIs, even in the same flow type. After extracting calibrated coefficients from the posterior distribution, they were computed in the deterministic solver and the results clearly showed what QoIs they were calibrated for. The results with calibrated coefficients show better agreement than GEKO-CC for both test cases.
The results obtained from the GEKO model using calibrated coefficients based on Bayesian inference show superior prediction for two types of compressible flow than the results from the default coefficients of the GEKO model.

Data Availability
Data are available on request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.   16 International Journal of Aerospace Engineering