Comparison of the Point-Collocation Non-Intrusive Polynomial (NIPC) and Non-Intrusive Spectral Projection (NISP) Methods for the γ − R θ Transition Model

Abstract

In the present work, a comparative study of two major non-intrusive polynomial chaos methods, Point-Collocation Non-Intrusive Polynomial Chaos (NIPC) and Non-Intrusive Spectral Projection (NISP), was conducted for the transitional $\gamma -R_{\theta }$ transitional model. Three multiple model coefficients, C_a2, C_e1, and C_e2 were considered with multiple random inputs with the assumption of uniform distributions with ±10% deviation. The target transitional flows were one around a flat plate and Aerospatiale A-airfoil. Deterministic solutions were obtained by employing the open source software OpenFOAM. The results of two methods were compared to the results of Monte Carlo simulation with 500 runs. The order convergence of the mean value and the standard deviation (STD) were compared in terms of the quantities of interest, drag and lift coefficients. Further, the most effective model coefficient for each transitional flow could be found through the calculation of the Sobol index.

1. Introduction

Nowadays, Computational Fluid Dynamics (CFD) has become an indispensable tool in the design and development of engineering products. The target is to reduce the time and cost for the new design of a prototype via computer simulations. Recently, the rapid progress of computation performance and reduction of computational cost have enabled the adoption of the Uncertainty Quantification technique (UQ) [1,2,3,4,5,6,7], which requires many simulation cases. One important part of UQ is uncertainty propagation, which is used to quantify the sensitivity and variation of the simulation results with consideration of single or multiple stochastic inputs such as the freestream velocity, model coefficients, and target geometry. A direct method for UQ is Monte Carlo simulation, but it is expensive to obtain converged solutions. Therefore, the Polynomial Chaos (PC) approach developed by Wiener [8,9] is proposed as an alternative method. Such methods can be divided into two types: intrusive and non-intrusive approaches. The intrusive methods require modification of the existing code into a Galerkin projection-based flow solver [10,11,12,13,14,15]. This extra work can be a big burden for CFD engineers. Whereas a non-intrusive method that enables the use of established code without any modification is considered as a promising approach in UQ and has been applied to various CFD fields by many researchers. Hosder and Walters [16,17] applied Point-Collocation Non-Intrusive Polynomial Chaos (NIPC) to selected CFD problems such as aerodynamic flows around the airfoil/wing or in a nozzle. They considered the uncertainty of the freestream velocity, eddy viscosity ratio, and model coefficients of the adopted turbulence models. They demonstrated that NIPC improved the accuracy of the polynomial chaos coefficient in NIPC approaches and reduced the computational expense by achieving the same level of the accuracy with a small number of sample points. Loeven [18] introduced a non-intrusive probabilistic collocation approach for efficient propagation of arbitrarily distributed parametric uncertainties. They also showed that the convergences depending on the order of approximation of the probabilistic collocation method and the Galerkin PC method are the same. In the UQ technique for turbulence model coefficients, Platteeuw et al. [19] applied the Probabilistic Collocation method to the model coefficients of the standard k-ε turbulence model. They found the necessary order for proper quantification of drag coefficients around a flat plate and NACA 0012 airfoil. Schaefer et al. [20] studied uncertainty quantification of three turbulence models: the Spalart–Allmaras model, Wilcox k-ω model, and Menter’s SST model for transonic wall-bounded flows. The significant model coefficients of three turbulence models were investigated through a sensitivity analysis.

The mechanism of transitional flow from laminar to turbulent flow is still a difficult problem in engineering and science fields. In particular, many researchers have developed various transition models based on Reynolds-Averaged Navier–Stokes (RANS) equations to apply to engineering problems to reduce the computational cost and time. Recently, Langtry and Menter [21,22,23,24] proposed a promising transition model based on local parameters of the intermittency γ and the momentum thickness Reynolds number Re_θ. This transition model has been applied to various engineering problems and the results such as the transition point and skin friction variation are superior to previous transition model results. However, in any transition or turbulence model, many assumptions and simplifications are required to close model coefficients with uncertainties. This factor was the motivation of the present work. We considered the uncertainties of the model coefficient of the $\gamma -R_{\theta }$ transition model and applied the UQ technique to the parameters.

In the present work, two non-intrusive methods, NIPC and the non-intrusive spectral projection method, were applied to the transitional flow around a flat plate and Aerospatiale A-Airfoil. The random inputs for UQ were the closure model coefficients of Menter’s $\gamma -R_{\theta }$ transition model and the quantities of the interests (QoIs) were the skin friction coefficient and drag/lift coefficients.

Even though the UQ methodologies used in the present work are classic ones in the UQ fields, we tried to access the application capability of two UQ methods with somewhat transitional algorithm in the transitional flow simulation; the first one was Non-Intrusive Spectral Projection (NISP) with the sampling of the quadrature rule and the second was point-collocation NIPC with the oversampling 2.

The theory of the non-intrusive method for UQ is described briefly in Section 2. In this section, we also explain the differences between the two methods, point-collocation NIPC and NISP. Next, the numerical method and Menter’s $\gamma -R_{\theta }$ transition model are presented in Section 3. In Section 4, the deterministic solutions for the flat plate and Aerospatiale A-Airfoil are validated with reference data and the stochastic results of two non-intrusive methods are compared to the Monte Carlo simulation with convergence with respect to the polynomial order and distribution of the probability density function of the outputs. Finally, Section 5 summarizes the main findings of the present study.

2. Non-Intrusive Method

This section introduces the theoretical background of the non-intrusive methods employed for the approximation of an output of a model involving random data, parameterized by a finite number of random variables $\xi (\theta )$ defined on a probability space $(\mathsf{\Theta },\Sigma ,P)$, and probability density $p_{\xi }$($\xi$). We denote $\Xi$ as the support of the density $p_{\xi }$ [15,25] and s represents the model output, the quantities of interest. The approximation of the mapping is obtained in terms of the basis of a random functional.

$$s:\xi \in \Xi \mapsto s\left(\xi \right)={\displaystyle \sum _i{s}_i{\mathsf{\Psi }}_i\left(\xi \right)}$$

(1)

The non-intrusive methods rely on a set of deterministic model resolutions, corresponding to some specific values or realizations of $\xi$ to construct the approximation s($\xi$).

The spectral modes (${\alpha }_k$) of the random variable ${\alpha }^{*}$ can be obtained by solving the linear system of equations given above. The mean (${\mu }_{{\alpha }^{*}}$) and the variance (${\delta }_{{\alpha }^{*}}^2$) of the solution can be obtained, as shown below.

$${\mu }_{{\alpha }^{*}}={\overline{\alpha }}^{*}\left(x\right)={\displaystyle {\int }_{\Xi }{\alpha }^{*}\left(x,\xi \right)p\left(\xi \right)d\xi }={\alpha }_0\left(x\right)$$

(2)

$$\begin{array}{ll}{\sigma }_{{\alpha }^{*}}^2& =\mathrm{var}\left[{\alpha }^{*}\left(x,\xi \right)\right]={\displaystyle {\int }_{\Xi }{\left({\alpha }^{*}\left(x,\xi \right)-{\overline{\alpha }}_0^{*}\left(x\right)\right)}^{2}p\left(\xi \right)d\xi }\\ & ={\displaystyle \sum _{i=1}^P\left[{\alpha }_{\mathrm{j}}^2\left(x\right)〈{\mathsf{\Psi }}_j^2\left(\xi \right)〉\right]}. \end{array}$$

(3)

2.1. Non-Intrusive Spectral Projection

We consider a model where both the input and output are real valued scalar random variables on a probability space $(\mathsf{\Theta },\Sigma ,P)$.

$$P_{\xi }(\xi )={\displaystyle \prod _{i=1}^N{p}_i\left({\xi }_i\right)}$$

(4)

Non-Intrusive Spectral Projection (NISP) [15] computes the projection coefficient of the random model output s($\xi$) on a finite dimensional stochastic subspace ${𝓈}^p$ of L₂$\left(\Xi ,p_{\xi }\right)$. The space equipped with the inner product $〈,〉$,

$$〈u,v〉={\displaystyle {\int }_{\Xi }u(\xi )v(\xi )p_{\xi }(\xi )d\xi },$$

(5)

is a Hilbert space. We denote (p + 1) as the dimension of ${𝓈}^p$ and let $\left\{{\mathsf{\Psi }}_0\left(\xi \right),\dots ,{\mathsf{\Psi }}_p\left(\xi \right)\right\}$ be an orthogonal basis of ${𝓈}^p$.

$$〈{\mathsf{\Psi }}_i,{\mathsf{\Psi }}_j〉={\displaystyle \int \dots {\displaystyle \int {\mathsf{\Psi }}_i\left(\xi \right){\mathsf{\Psi }}_j\left(\xi \right)p_{\xi }\left(\xi \right)d{\xi }_1\dots d{\xi }_N=〈{\mathsf{\Psi }}_i,{\mathsf{\Psi }}_j〉{\delta }_{ij}}}.$$

(6)

Non-Intrusive Spectral Projection provides a decomposition of (s)^|| into a finite generalized Polynomial Chaos (gPC) basis $\left\{{\mathsf{\Psi }}_0\left(\xi \right),\dots ,{\mathsf{\Psi }}_p\left(\xi \right)\right\}$. Thus, it is the truncated Polynomial Chaos Expansion (PCE) of the probability space. The spectral projection of the output is shown below.

$$s^{\left|\right|}\left(\xi \right)={\displaystyle \sum _{k=0}^P{s}_k{\mathsf{\Psi }}_k(\xi )\in }{𝓈}^p$$

(7)

Since $s^{\left|\right|}\left(\xi \right)$ is the projection of the output onto ${𝓈}^p$, for each index $k\in \left\{0,\dots ,n\right\}$, the coefficients are given by

$$s_k=\frac{〈s,{\mathsf{\Psi }}_k〉}{〈{\mathsf{\Psi }}_k,{\mathsf{\Psi }}_k〉}, k = 0, 1, \dots ,$$

(8)

such that the following is satisfied.

$$s-s^{\left|\right|}\perp {𝓈}^p\iff 〈s-s^{\left|\right|},{\mathsf{\Psi }}_k〉=0 \mathrm{for} 0\le k\le P$$

(9)

When the generalized Polynomial Chaos (gPC) basis is selected, the values $〈{\mathsf{\Psi }}_k,{\mathsf{\Psi }}_k〉$ are known analytically and, therefore, only the computation of the scalar product remains:

$$〈{\mathsf{\Psi }}_k,{\mathsf{\Psi }}_k〉s_k=〈s,{\mathsf{\Psi }}_k〉={\displaystyle {\int }_{\Xi }s({\xi }_1,\dots ,{\xi }_N){\mathsf{\Psi }}_k{p}_1({\xi }_1)}\dots p_N({\xi }_N)d{\xi }_1\dots d{\xi }_N.$$

(10)

Here, $\Xi \subset R$ is the support of the probability density function $p_{\xi }$ of the basic random variable $\xi$. The multidimensional integral of a function $f(\xi )\equiv s(\xi ){\mathsf{\Psi }}_k(\xi )$ over $\Xi$ with a non-negative weight $p_{\xi }$ has a product form, which can be shortened as shown below.

$$\mathsf{\gamma }f={\displaystyle \underset{\Xi }{\int }f(\xi )p_{\xi }(\xi )d\xi }$$

(11)

When the dimensionality N of the integration is not too large, as is the case for many problems in practical engineering fields, deterministic cubature can be effectively used. A cubature is an approximation of the multidimensional integral γ$f$ as a discrete sum:

$$\mathsf{\gamma }f\approx {\displaystyle \sum _{i=1}^{N_Q}f({\xi }^{(i)})w^{(i)}},$$

(12)

where ${\xi }^{(i)}\in \Xi$ and $w^{(i)}\in R$, i = 1, …, N_Q, are the nodes and weights of the chosen quadrature formula. Various quadrature rules can be used by means of a non-negative weight function $p_{\xi }$. N_Q is the number of model resolutions.

In further sections, the set ${\left\{{\xi }^{(i)}\right\}}_{i=1}^{N+1}$ is called a sampling of $\xi$ due to the very definition of the integral in Equation (12) and the quadrature node can be understood as particular realizations of a basic random variable $\xi$.

By combining Equations (8), (10), and (12), it becomes clear how NISP relies on a set of deterministic model resolutions, corresponding to some specific realizations of $\xi$. Along this line, a deterministic simulation code can be used as a black-box, which associates each realization of the parameters to the corresponding model outputs.

2.2. Point-Collocation Non-Intrusive Polynomial Approach

Unlike the Non-Intrusive Spectral Projection (NISP) method discussed above, this method does not aim to determine the projection of the model output s($\xi$) on a predefined stochastic subspace, but instead relies on the interpolation based on Equation (7). Then, P + 1 vectors $\left(\overrightarrow{\xi }={\left\{{\xi }_1,{\xi }_2,\dots ,{\xi }_n\right\}}_k,k=0,1,2,\dots ,P\right)$ are chosen in the random space for a given PC expansion with P + 1 modes and the deterministic simulation code is evaluated at this point, which has a specific probability distribution. The discrete sum is taken over the number of output modes [15],

$$P+1=\frac{(n+p)!}{n!p!},$$

(13)

where p and n are the order of polynomial chaos and the number of the random input dimension, respectively. The PC model obtained via Equation (13) is formed by polynomials up to a maximum degree P (also called order of the PC model): this truncation strategy is referred to as a “total order expansion”. A linear system of equations can be obtained [16,17,18,25,26,27], as shown below.

$$\left[\begin{array}{cccc}{\mathsf{\Psi }}_0\left(\overrightarrow{{\xi }_0}\right)& {\mathsf{\Psi }}_1\left(\overrightarrow{{\xi }_0}\right)& \cdots & {\mathsf{\Psi }}_p\left(\overrightarrow{{\xi }_0}\right)\\ {\mathsf{\Psi }}_0\left(\overrightarrow{{\xi }_1}\right)& {\mathsf{\Psi }}_1\left(\overrightarrow{{\xi }_1}\right)& \cdots & {\mathsf{\Psi }}_p\left(\overrightarrow{{\xi }_1}\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\mathsf{\Psi }}_0\left(\overrightarrow{{\xi }_p}\right)& {\mathsf{\Psi }}_1\left(\overrightarrow{{\xi }_p}\right)& \cdots & {\mathsf{\Psi }}_p\left(\overrightarrow{{\xi }_p}\right)\end{array}\right]\left\{\begin{array}{c}{\alpha }_0\\ {\alpha }_1\\ ⋮\\ {\alpha }_p\end{array}\right\}=\left\{\begin{array}{c}{\alpha }^{*}\left(\overrightarrow{{\xi }_0}\right)\\ {\alpha }^{*}\left(\overrightarrow{{\xi }_0}\right)\\ ⋮\\ {\alpha }^{*}\left(\overrightarrow{{\xi }_0}\right)\end{array}\right\}$$

(14)

We solved the linear system of Equations (14), which requires P + 1 deterministic function evaluations. If more than P + 1 samples are chosen, then the over-determined system of equations should be solved using the Least Squares method. Hosder et al. [16] investigated the effects on the results by the number of collocation points in a systematic way through the introduction of a parameter, the oversampling rate n_p defined below.

$$n_p=\frac{number of samples}{P+1}$$

(15)

They used oversampling rates of n_p = 1, 2, 3, and 4 for the solution of a stochastic model problem with multiple uncertain variables to study the effect of the number of collocation points on the given accuracy of the polynomial chaos expansions. They suggested that a number of collocation points corresponding to two times the minimum number required, which means that an oversampling rate n_p of 2 yields a better approximation to the statistic at each polynomial degree.

The open source packages Chaospy [28] and Scilab [15,29,30] were used for UQ analysis in the present work. The Chaospy [28] software package provides various stochastic tools for UQ analysis for the point-collocation NIPC method. Scilab software [15] was developed for spectral projection based on the Gaussian quadrature.

3. Numerical Methods and Menter’s $\gamma -R_{\theta }$ Transition Model Formulation

In the present work, the open source code OpenFOAM V.5.0 [31] was adopted for the deterministic flow solver. The steady incompressible Navier–Stokes governing equations were discretized spatially at the second order and the SIMPLE algorithm was adopted for continuity.

Menter proposed the $\gamma -R{e}_{\theta }$ transitional model [21] based on the fully turbulent SST k-$\omega$ model. It includes two transport equations for intermittency, $\gamma$, and for the transported transition momentum thickness Reynolds number ${\widetilde{\mathrm{Re}}}_{{}_{\theta t}}$. The dimensionless intermittency and Reynolds number based on the momentum thickness are formulated as follows.

$$\rho \frac{D\gamma }{Dt}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_t}\right)\frac{\partial \gamma }{\partial x_j}\right]+F_{length}{c}_{a1}\rho S{\left(\gamma F_{onset}\right)}^{0.5}(1-c_{e1}\gamma )+c_{a2}\rho \mathsf{\Omega }\gamma F_{turb}(1-c_{e2\gamma })$$

(16)

$$\rho \frac{D\tilde{\mathrm{R}}{e}_{\theta t}}{Dt}=\frac{\partial }{\partial x_j}\left[{\sigma }_{\theta t}(\mu +{\mu }_t)\frac{\partial \widetilde{{\mathrm{Re}}_{\theta t}}}{\partial x_j}\right]+c_{\theta t}\frac{{(\rho U)}^2}{500\mu }({\mathrm{Re}}_{\theta t}-D\tilde{\mathrm{R}}{e}_{\theta t})(1-F_{\theta t})$$

(17)

The modeled transport equations are controlled by the following functions:

$$F_{onset}=\max (F_{onset2}-F_{onset3},0),F_{turb}=\exp [-{(R_T/4)}^4],$$

(18)

$$F_{onset2}=\min [\max (F_{onset1},F_{onset1}^4),2],F_{onset1}=\frac{{\mathrm{Re}}_v}{2.193{\mathrm{Re}}_{\theta c}},$$

(19)

$$F_{onset3}=\max (1-{(R_T/2.5)}^3,0),$$

(20)

$$F_{\theta t}=\min \left\{\max \left[F_{wake}.\exp \left[-{\left(\frac{U^2}{375\mathsf{\Omega }v\tilde{R}{e}_{\theta t}}\right)}^4\right],1-{\left(\frac{c_{e2}\gamma -1}{c_{e2}-1}\right)}^2\right],1\right\},F_{wake}=\exp \left[-{\left(\frac{{\mathrm{Re}}_{\omega }}{{10}^5}\right)}^2\right],$$

(21)

with the following parameters:

$${\mathrm{Re}}_v=\frac{\rho S{y}^2}{\mu },R_T=\frac{\rho k}{\mu \omega },{\mathrm{Re}}_{\omega }=\frac{\rho \omega y^2}{\mu },{\mathrm{Re}}_v=\frac{\rho S{y}^2}{\mu },R_T=\frac{\rho k}{\mu \omega },{\mathrm{Re}}_{\omega }=\frac{\rho \omega y^2}{\mu }$$

(22)

The details of the intermittency equation are presented in a previous work [21]. Both F_length and $R_{\theta c}$ in the above equation are empirical correlations and functions of the transition momentum-thickness Reynolds number $R_{\theta t}$, which is also a function of the turbulence intensity $Tu$ and the streamwise pressure gradient ${\lambda }_{\theta }$ in the free flow. The model constants are c_a1 = 2.0, c_e1 = 1.0, c_a2 = 0.06, c_e2 = 50.0, $c_{\alpha }$ = 0.5, ${\sigma }_{\gamma }$ = 1.0, ${\sigma }_{\theta t}$ = 2.0 and $c_{\theta t}$ = 0.03.

For predicting separation induced transition, the following modification is utilized.

$${\gamma }_{sep}=\min (2.0\max (0,\frac{{\mathrm{Re}}_v}{3.235R_{\theta c}}-1)F_{reattach},2)F_{\theta t}$$

(23)

$$F_{reattach}=e^{-{\left(\frac{R_T}{20}\right)}^4}$$

(24)

The effective value of $\gamma$ is thus obtained as follows.

$${\gamma }_{effective}=\max (\gamma ,{\gamma }_{sep})$$

(25)

The empirical correlation used in this model is based on the pressure gradient parameter and turbulent intensity defined as follows.

$${\lambda }_0=\frac{{\theta }^2}{v}\cdot \frac{dU}{ds},K=\frac{v}{U^2}\cdot \frac{dU}{ds},Tu=\frac{100{(2k/3)}^{1/2}}{U},\frac{dU}{ds}=\frac{\partial u_j}{\partial x_j}\frac{u_i{u}_j}{U^2}$$

(26)

${\mathrm{Re}}_{\theta }$ is defined as follows.

$${\mathrm{Re}}_{\theta t}=\frac{\rho {\theta }_t{U}_0}{\mu }$$

(27)

The following correlation equations were defined by Langtry [21].

$${\mathrm{Re}}_{\theta t}=\left\{\begin{array}{l}\left[1173.51-589.428Tu+\frac{0.2196}{T{u}^2}\right]F(\lambda ,Tu);\hspace{1em}Tu\le 1.3\\ 331.5{[Tu-0.5658]}^{-0.671}F(\lambda ,Tu);\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}Tu>1.3\end{array}\right.$$

(28)

$$F(\lambda ,Tu)=\left\{\begin{array}{l}1+e^{-{(2Tu/3)}^{1.5}\cdot (12.986\lambda +123.66{\lambda }^2+405.689{\lambda }^3);\hspace{1em}\hspace{1em}\lambda \le 0}\\ 1+0.275e^{-2Tu}\cdot (1-e^{(-35\lambda )});\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\lambda >0\end{array}\right.$$

(29)

4. Problem Description and Results

In the $\gamma -R_{\theta }$ transition model, two transport equations with an extra nine model coefficients should be solved additionally based on the SST k-$\omega$ turbulent model. If a total of nine model coefficients are considered for multiple random inputs, tremendous computational time and cost are required. For example, when the order of the polynomial chaos p is 2 or 3, the required number of random inputs is 19,683 or 262,144, respectively in the Non-Intrusive Spectral Projection (NISP) method. First, the sensitivity analysis of each model coefficient was conducted with ±10% deviation for the flat plate and the three most effective model coefficients (C_a2, C_e1, and C_e2) were obtained for the multiple random inputs. The deviation and distribution of the model coefficients were set based on References [16,17,18].

Table 1. The sensitivity analysis of each model coefficient.

Coefficient	CD of Lower Boundary (−10%)	CD of Upper Boundary (10%)	Difference Percent (%)
Ca2	0.00391	0.00410	4.60
Ce1	0.00373	0.00421	11.4
Ce2	0.00397	0.00406	2.22
Ca1	0.00395	0.00402	1.74
Cα	0.00401	0.00405	0.99
${\delta }_{\gamma }$	0.00399	0.00401	0.49
${\sigma }_{\theta t}$	0.00389	0.00392	0.76
$C_{\theta t}$	0.00402	0.00398	1.0

presents the results of the sensitivity analysis which show the difference with the reference value of the deterministic solution. To carry this out, each model coefficient was assumed to have a uniform distribution with the fixed constant from the original model [21] and ±10% deviation. Then, we applied the UQ technique with the single random input to the target problem. We calculated the variance depending on each random model coefficient and determined the model coefficients that show the largest variance of the quantities of interest. Considering the simulation time and cost, the three top model coefficients c_a2, c_e1, and c_e2 were selected and applied to the NIPC and NISP methods with multiple random inputs. The following

Table 2. The uniform distributions of the three model coefficients ±10% deviation.

Coefficient	Original Value	Lower Boundary (−10%)	Upper Boundary (10%)
Ca2	0.06	0.054	0.066
Ce1	1.0	0.9	1.1
Ce2	50.0	45.0	55.0

shows the probability distributions of the three coefficients.

From the probability distributions of the three coefficients shown in

Table 3. The total numbers of selected points for Monte Carlo, NIPC, and NISP for three chosen points (n is the dimension and p is the order of stochastic computation).

n = 3	Monte Carlo MC	Point-Collocation NIPC	Spectral Project NISP
Order	-	$2\frac{(n+p)!}{n!p!}$	${(p+1)}^n$
2		20	27
3	500	40	64
4		70	125
5		112	216

, in our calculations, we used multi-dimensional Legendre polynomials, which are orthogonal in the interval for each random dimension. Stochastic solutions to the target problems were obtained using three approaches: Monte Carlo with 500 samples, NIPC, and NISP methods. For NIPC, the oversampling rate $n_p$ was set to 2, as recommended by Hosder et al. [16]. The Gaussian quadrature rule was adopted to select the sample points by the specified order, p + 1 for NISP.

The $\gamma -R_{\theta }$ transition model introduced by Menter [21] was applied to two target problems: transitional flow over a flat plate and flow around Aerospatiale A-airfoil [32,33,34]. These are well-known validation problems for validation of CFD algorithms and transition/turbulence models. Experimental data are presented on the site of turbulence modeling resources by NASA [34]. The deterministic solutions of two transitional flows were validated by comparison with reference data.

4.1. Transitional Flow over a Flat Plate

The first case was the transitional flow over a flat plate and the deterministic solution of this case was validated using the reference data of Langtry and Menter [21]. Figure 1 shows the computational mesh for the flat plate and the boundary conditions of the computational domain over the flat plate. The freestream turbulence intensity and eddy viscosity ratio values were set to $T_u=3.3\%$ and ${\mu }_t/\mu =10$ (Savill 1996) [33], respectively. A mesh with 43,000 cells was adopted and the mesh was clustered to the wall for $y^{+}\le 1$ in

Table 4. Information regarding the inflow conditions and grid system for flow over a flat plate.

Inflow Condition	$U_{\infty }$	5.4 m/s
	$T_u$	3.3%
	${\mu }_t/\mu$	10.0
Grid System	(Lx,Ly)	(3.2c,1.0c)
	(Nx,Ny)	(230,100)
	${\mathrm{Re}}_L$	1 × 106 (L = 0.3 m)

. The corresponding Reynolds number based on the free stream velocity (5.4 m/s) and length of the flat plate (0.3 m) was 1.0 × 10⁶. Figure 2 shows the comparison of the skin friction $\left(C_f\right)$ profiles predicted by the deterministic solution with the $\gamma -R_{\theta }$ transition model and experimental data [33,34]. The skin friction coefficient determined by the present deterministic solution could predict the transition point and developing turbulent flow well. The slight early recovery from laminar to turbulent flow and underestimation of the skin friction right after transition agreed well with previous simulation results [21].

In this study, the drag coefficient was considered as the quantity of interest of UQ. Before discussing the stochastic results of UQ, the predicted drag coefficient obtained by the deterministic solution with the originally fixed model coefficients was 0.00401.

Table 5. The mean and standard deviation of the drag coefficient values for flat plate.

Order	Mean_NIPC	STD_NIPC	Mean_NISP	STD_NISP	Mean_MC	STD_MC
2	0.00413298	0.0001857	0.00411269	0.0001878	0.004128139	0.00017982
3	0.00412933	0.0001811	0.00412791	0.0001769
4	0.00412997	0.0001826	0.00412113	0.0001819
5	0.00412992	0.0001815	0.00412758	0.0001775

shows the drag coefficients predicted by the uncertainty quantification technique with the assumptions of uniform distribution of three model coefficients, c_a2, c_e1, and c_e2, and three methods: two methods are non-intrusive polynomial chaos methods (point-collocation NIPC and non-intrusive spectral [rojection) and a Monte Carlo method with 500 simulations. The mean value and standard deviation (STD) of the drag coefficient of all three methods are presented with respect to the polynomial order.

The mean drag coefficients approached approximately 0.00413, which showed 2.9% difference compared with the deterministic solution, 0.00401. The standard deviation, which could not be predicted by the deterministic simulation, was approximately 0.00018, which was 4.4% of the mean value.

The percent errors of the mean and STD values were, respectively, 0.043% and 0.9% for point-collocation NIPC and 0.014% and 1.29% for NISP when the two methods with fifth-order were compared to those in the Monte Carlo method. Figure 3 shows the convergence of the mean and STD of the drag coefficient according to the polynomial order of the two methods. The long dashed line is the reference value, which corresponds to the Monte Carlo results. The point-collocation NIPC showed faster convergence of the mean and STD than NISP and consistent values except for the case with second order. However, in NISP, the results of the odd number orders (third and fifth orders) are close to the Monte Carlo results and the results of the even number orders (second and fourth orders) have some discrepancies compared to the reference values.

Figure 4 shows the detailed probability density function of the drag coefficient, which was generated from the calculated polynomial coefficients based on 10,000 sampling points by the Latin-Hypercube sampling method. In all cases, the peaks of the probability density function were biased toward lower values with the given interval but the value of the peak showed small differences depending on the polynomial order. As the order increased, the shape of this function had a similar pattern but there was a small difference towards higher drag coefficient values.

4.2. Transitional Flow over Aerospatiale-A Airfoil

In this section, the transitional flow around the Aerospatiale A-airfoil is presented under the conditions of the maximum lift angle of attack ($\alpha ={13.1}^{\circ }$) and ${\mathrm{Re}}_c=2.07\times {10}^6$ based on the freestream velocity and chord length as shown in Figure 5. The information of the geometry and the experimental results for the Aerospatiale-A airfoil is available from the Office National d’Etudes et Recherches Ae1rospatiales F1 (ONERA F1) wind tunnel (Chaput 1997) [32]. The inflow conditions and grid system are summarized in

Table 6. Information about the inflow conditions and grid system for flow over Aerospatiale A-airfoil.

Inflow Condition	$U_{\infty }$	32.8 m/s
	$T_u$	0.2%
	${\mu }_t/\mu$	10.0
	$\alpha$	${13.1}^{\circ }$
Grid System	(Lx,Ly)	(45c,30c)
	(Nx,Ny)	(350,120)
	${\mathrm{Re}}_c$	$2.07\times {10}^6$

. The turbulent intensity and eddy viscosity ratio were set to $T_u=0.2\%$ and ${\mu }_t/\mu =10$ [20], respectively. The length of the upstream part and downstream part were 20c and 25c and the height of the computational domain was 30c. At the upstream and side parts, the velocity inlet boundary condition was given and, at the downstream, the Neumann boundary condition Was set for the outlet. To satisfy $y^{+}\le 1$ for correct resolving of the turbulent boundary layer, the distance of the first mesh off the wall was ${10}^{-5}$ c and the growth ratio was 1.2.

Figure 6 shows the skin friction coefficient distributions over the upper wall of the airfoil predicted by the deterministic solver with the $\gamma -R_{\theta }$ transition model as well as the comparison with experimental data [32]. The transition point was predicted to be between 10% and 15% of the chord length, which agreed well with other simulation results [35]. The behavior with developing turbulent flow after transition showed consistent results with experimental data.

Further, the coefficients of the drag and lift forces were calculated and compared to experimental data, F1 and F2 [32] as show in

Table 7. Lift and drag at $\alpha ={13.1}^{\circ }$

Model	${\mathit{C}}_{\mathit{L}}$	${\mathit{C}}_{\mathit{D}}$
$\gamma -R_{\theta }$	1.390	0.0206
F1 exp	1.562	0.0208
F2 exp	1.515	0.0308

. Unfortunately, there were large discrepancies between the experimental data, F1 and F2, especially for the drag coefficient. When the measurement sensitivity of this transitional flow was considered, the present simulation showed reasonable results similar to previous simulation results [35,36]. It was confirmed that the present simulation adopts adequate methodologies and grid system.

The mean value and STD of the drag and lift coefficients of the Aerospatiale A-Airfoil were calculated using the point-collocation NIPC, NISP, and Monte Carlo methods. The number of samples of Monte Carlo was 500, which was the same as in the previous case.

Table 8. The mean and standard deviation of the drag coefficient (C_D) values for A-airfoil.

Order	Mean_NIPC	STD_NIPC	Mean_NISP	STD_NISP	Mean_MC	STD_MC
2	0.020561	0.0001252	0.020598	0.0000895	0.02057489	0.0001119
3	0.020570	0.000111	0.020508	0.0001090
4	0.020551	0.000113	0.020587	0.0001060
5	0.020569	0.000112	0.020581	0.0001130

shows the drag coefficient with respect to the polynomial order.

When the polynomial order is 5, the mean drag coefficient was predicted to be 0.020569 by PC-NIPC and 0.020581 by NISP. When these were compared with the deterministic solution of the original model coefficients, it was found that the drag coefficient was not as biased as in the flat plate case.

The percent errors of the mean and STD values were 0.029% and 0.09% for point-collocation NIPC and 0.03% and 0.97% for NISP when the two methods with fifth-order were compared to the Monte Carlo method. When the convergence of the drag coefficient with respect to the polynomial order was checked, the overall pattern was similar to the previous case for the drag coefficient of the flat plate flow. The mean value of the drag coefficient already converged within 0.5% (ΔC_D ~ 0.0001) at the second order, whereas the STD converged after the second polynomial order as shown in Figure 7. As the polynomial order increased, the mean and STD approached the values obtained by the Monte Carlo method (long dashed line).

The distributions of the probability density functions of the drag coefficients are plotted with variation of the order in Figure 8. In this case, the distributions of the PDF revealed different patterns depending the adopted method and the order of the polynomial. In the case of point-collocation NIPC, the distribution of the case with second order was very different form other order cases and had similar patterns starting from the third order. NISP had a similar trend as point-collocation and the results were similar starting from the fourth-order case. Interestingly, there were two peaks in the distribution of the PDF in high-order cases for the two methods. The transitional flow simulation of Aerospatiale A-Airfoil was more sensitive to used grid, convergence criteria and the adopted scheme than the flat plate simulation. This means that the quantities of interest, such as drag and lift coefficients, were sensitive to the computational conditions (grid, convergence criteria and so on) and sampling points that were selected at the given order of polynomial. We presume that the sampling points were insufficient to resolve the exact probabilistic distribution in the lower Orders 2 and 3.

Table 9. The mean and standard deviation values of the lift coefficient (C_L) for A-airfoil.

Order	Mean_NIPC	STD_NIPC	Mean_NISP	STD_NISP	Mean_MC	STD_MC
2	1.39999	0.017730	1.397522	0.017321	1.399155	0.016661
3	1.39955	0.016821	1.399802	0.016343
4	1.39949	0.016969	1.398551	0.017026
5	1.39967	0.016804	1.399507	0.016583

shows the mean and STD of the lift coefficient obtained by the three methods. The mean values of all methods were predicted to be approximately 1.399, whereas the deterministic solution was 1.390.

Similar to the convergence of the drag coefficient, the mean value of the lift coefficient converged early from the second order but the STD showed a large difference when the order was 2. At the fifth order, there was a 0.85% difference in point-collocation NIPC and 0.46% in NISP as shown in Figure 9. Even though the convergence rate of point-collocation NIPC was faster than that of NISP, the accuracy at the highest order (fifth order) was higher in NISP than point-collocation NIPC. The difference between predicted and the deterministic solution value resulted because the distribution of PDF was biased toward the positive difference value, as confirmed in Figure 10.

In the lift coefficient, similar patterns of the PDF were obtained depending on the method and the order. Even though the peak of PDF was located at a lower value, the flat distributions of the lift coefficients ranged to the value of 1.43 or 1.44 in every case. This broad distribution towards a higher lift coefficient made the predicted mean lift coefficient move to a higher value than the deterministic solution.

4.3. Comparison of the Two Methods

Figure 11 describes the distribution of the PDF of point-collocation NIPC and NISP for the fifth-order polynomial and Monte Carlo method with 500 sampling points. In the case of flat plate, the peak values of the drag coefficients, 0.0041, agreed well in the three methods, NISP, PC-NIPC and Monte Carlo. This peak was located at a higher value than that (0.00401) of the deterministic solution. In the drag coefficients of Aerospatiale A-Airfoil simulation, there were two peaks in the PDF and the left peak was lower than the right one in every method. The predicted peak values in NISP and PC-NIPC methods were approximately 0.0206 with 0.03% error with that of Monte Carlo method. The pattern of the distribution of the lift coefficient in A-Airfoil simulation was similar to that of the drag coefficient in the flat plate simulation. The peak value of NISP showed smaller difference with 0.46% than that of PC-NIPC (0.86%) when compared with the Monte Carlo result. Even though the number of samplings (500) in the present work was not sufficient to get the converged solution of MC, the overall pattern distribution and the number of the peak were consistent between MC and two adopted methods. When the number of sample points for the UQ technique was considered, even though point-collocation NIPC required fewer samples than NISP (Table 3), similar accuracy and PDF as the mean and STD of the quantities of interest could be obtained. In particular, the number of random variables increased and the point-collocation NIPC was useful from a computational cost and time point of view.

Sensitivity analysis of three model coefficients, c_a2, c_e1, and c_e2, was carried out using the Sobol index of the two methods (

Table 10. Sobol index of three model coefficients for point-collocation NIPC.

Sobol Index	Ca2	Ce1	Ce2
CD Flat plate	0.01934	0.983	0.01955
CD A-airfoil	0.32256	0.8963	0.15859
CL A-airfoil	0.05347	0.9793	0.03839

and

Table 11. Sobol index of three model coefficients for NISP.

Sobol Index	Ca2	Ce1	Ce2
CD Flat plate	0.0034	0.9945	0.0022
CD A-airfoil	0.0287	0.9554	0.0281
CL A-airfoil	0.00024	0.9999	0.000234

). The Sobol indices provided quantitative information regarding how much effect the corresponding random parameters had on quantities of interest. In the point-collocation NIPC and NISP methods, the most effective model coefficient was calculated as C_e1, whose Sobol index was generally over 90% in all outputs of both the flat plate and A-airfoil cases. The model coefficient C_e1 was included in the source term of the transport equation of the intermittency and was related to the F_legnth and F_onset terms. The other two variables, C_e2 and C_a2, showed similar effects on the results of less than 5% of Sobol indices. This means that more careful calibration of the model coefficient Ce₁ was required compared to the other coefficients for better prediction of transitional flow.

Figure 12 displays the distribution of the skin friction coefficient from point-collocation NIPC at the fifth order. It was confirmed that the change of the model coefficients had a significant effect on the variation of the skin friction coefficient at the location of the flat plate. There was a small STD before the transition point (i.e., laminar flow region) and then a large STD after the transition (i.e., developing turbulent flow region). This means that the model closure coefficients had a large effect on the turbulent flow region of the deterministic solution.

5. Conclusions

In this study, we compared point-collocation NIPC and NISP methods for the $\mathsf{\gamma }-{\mathrm{Re}}_{\mathsf{\theta }}$ transitional model [21]. Three model coefficients, Ca₂, Ce₁, and Ce₂, were considered as multiple random inputs with the assumption of a uniform distribution with ±10% deviation. The quantities of interest were the drag and lift coefficients in transitional flow around a flat plate and Aerospatiale A-Airfoil.

The point-collocation NIPC had the advantage of having fewer samples for uncertainty quantification analysis, even though the oversampling ratio was set to 2, followed by the NISP method. It was confirmed that the point-collocation NIPC is useful when multiple random inputs are considered. Convergence with respect to the polynomial order was faster in point-collocation NIPC than NISP. In addition, the mean value of the output was obtained with a lower order than STD. This means that a higher polynomial order was necessary to obtain a converged STD. At the highest polynomial order of 5 evaluated in this study, NISP demonstrated slightly more accurate results when compared to the Monte Carlo simulation. However, the accuracy of NISP depended on whether it was an even or odd polynomial order. Through analysis of the Sobol index, the most effective model coefficient of the $\gamma -R{e}_{\theta }$ transitional model was Ce₁, which is related to F_length and F_oneset in the intermittency equation, by approximately over 90%.

In summary, when multiple random inputs were considered, point-collocation NIPC could provide reasonable stochastic values such as the mean and STD with fewer samples than NISP. If only one or two variables were regarded as random inputs, NISP seemed to be better at a high order.