Uncertainty Quantification for Multidimensional Correlated Flow Field Responses

With the ongoing advancements in mathematical modeling, numerical algorithms, grid techniques, and high-performance computing, computational fluid dynamics (CFD) has emerged as an increasingly critical component in numerous engineering disciplines, including aerospace, energy and power generation, transportation, and beyond. This growth has been accompanied by a corresponding increase in attention to the credibility of CFD. Nevertheless, uncertainties inherent to CFD persist, encompassing models, parameters, numerical solutions, and more [1], prompting concerns regarding its credibility. Therefore, a rigorous and systematic quantification of uncertainties within CFD is of utmost importance in assessing and enhancing its credibility.

Parameters remain crucial uncertain factors in CFD. Owing to cognitive constraints or the intrinsic randomness, model parameters or the inflow environment can exhibit a certain degree of uncertainty. To address this issue, scholars have devised a range of methodologies, such as Monte Carlo (MC) simulations, polynomial chaos [2] and evidence theory [3,4], to quantify the uncertainty introduced by parameters, achieving notable success in engineering applications. In recent years, scholars have been exploring the integration of machine learning algorithms into the field of uncertainty quantification (UQ). By creating a low-cost surrogate model for the response, it can be utilized as a substitute for the original simulator, enabling the utilization of Monte Carlo techniques for UQ tasks. Well-regarded surrogate models in the literature include Gaussian processes [5,6], radial basis functions [7], and deep neural networks [8], among others.

Nevertheless, most of these techniques have been limited to single-output scenarios. In the context of CFD, the desired outputs encompass not only individual variables but also time- or space-varying flow field variables, such as wall pressure coefficient distribution and unsteady aerodynamic forces of aircraft. The potential correlation between flow field variables at different locations or times results in high-dimensional outputs, even reaching tens of thousands. Currently, three primary approaches exist for modeling multidimensional correlated responses: (1) the multiresponse model [9,10], which treats the output of interest as a vector and utilizes a multi-output Gaussian model for modeling. This approach necessitates estimating the correlation between output dimensions, leading to a substantial increase in the number of model parameters to be estimated. (2) The joint space–time model [11–13], which equates time/space variables to model input variables, converting them into scalar outputs. This method results in a training sample size proportional to the product of the original sample size and the output dimension, potentially leading to excessive training sample size. (3) The individual surrogate model [14], which models each output dimension independently, disregarding the correlation among them. This approach does not guarantee modeling accuracy and may incur substantial modeling costs. Therefore, there is an urgent need to develop efficient and reliable uncertainty propagation modeling techniques for multidimensional correlated responses.

In recent years, flow field reduction techniques exemplified by proper orthogonal decomposition (POD) have gained widespread application in the analysis of flow characteristics. Scholars have employed POD to investigate turbulence flow features [15], aircraft vortex structures [16], and the evolution of coherent structures in separated flows [17]. The obtained results attest to the robust representational capabilities of POD in extracting unsteady flow field characteristics, enabling accurate reconstruction of the original flow field using a limited number of basis functions. In the analysis of unsteady flow features, time serves as a variable parameter, whereas in the quantification of flow field uncertainty, model parameters assume the role of variable parameters, sharing conceptual similarities. This provides a novel avenue for quantifying flow field uncertainty under the influence of uncertain parameters.

This paper presents a modeling approach for uncertainty propagation based on POD and artificial neural networks (ANN), specifically targeting multidimensional correlated responses in CFD. By applying POD, the multidimensional correlated responses are projected onto an orthogonal basis function space, enabling a bidirectional representation between the responses and the basis function vectors. By utilizing ANN, a predictive model is constructed between model uncertainty parameters and basis function vectors, enabling the prediction of multidimensional correlated responses under arbitrary uncertain model parameters via the bidirectional representation. To effectively identify model parameters with substantial impacts on multidimensional correlated responses, a multi-output global sensitivity analysis method grounded in covariance decomposition is adopted.

The paper is organized as follows: Section 2 introduces the proposed uncertainty propagation model and global sensitivity analysis method for multidimensional correlated responses. Section 3 presents the case studies, which are the low-speed flow around a NACA0012 airfoil and the transonic flow around a M6 wing. Section 4 contains the analysis results of the cases, including the prediction accuracy analysis of the entire method, statistical information of the flow field, and parameter sensitivity results. Finally, Sec. 5 discusses the conclusion and future directions for further research.

This paper employs a surrogate model to characterize the propagation of uncertainty. Given the multidimensional and correlated nature of the responses, it is essential to model the output holistically. Although techniques such as ANN are capable of modeling high-dimensional outputs, the substantial number of model parameters to be estimated can impact training time. Therefore, this article initially utilizes POD to reduce the dimensionality of the flow field. The flowchart of the whole approach is shown in Fig. 1. Initially, we employ Latin hypercube sampling to sample the uncertain model parameters within their defined parameter space, generating the input data for training samples. Next, CFD simulations are performed to generate the multidimensional flow field responses for each training sample. We utilize the POD technique to reduce the dimensionality of these responses, creating a reduced-dimension representation. In the fourth step, a surrogate model (here using ANN) is used to establish a prediction model between the uncertain model parameters and the reduced-dimension responses. This surrogate model can predict reduced-dimension responses for arbitrary inputs within the parameter space. Since the POD method is a bidirectional representation between the original multidimensional correlated responses and the reduced-dimension responses, the corresponding multidimensional correlated responses can be predicted. Based on the whole model, large-scale random sampling can be carried out to analyze the statistical information regarding response quantities, such as mean values, confidence intervals, and other metrics. Additionally, we can assess the global sensitivity of input parameters to identify which parameters have the greatest impact on overall response.

The first part of this section will introduce the key algorithms in the uncertainty propagation model, including POD and ANN algorithms, and the second part will introduce the global sensitivity analysis method.

Proper orthogonal decomposition is a powerful mathematical technique used for analyzing and reducing the complexity of high-dimensional systems. This method is widely applied in various fields such as fluid dynamics, structural dynamics, and signal processing, enabling efficient extraction of important characteristics, reduction of computational costs, and enhancement of simulation efficiency. POD provides a valuable tool for the analysis and optimization of complex systems.

Proper orthogonal decomposition method was proposed by Lumley [15], and its basic principle is to find a new set of basis for the space spanned by a group of vectors in high-dimensional space such that the projection of the original vectors onto the set of basis is as large as possible. It involves projecting data onto a set of orthogonal basis functions, allowing for the extraction of dominant modes and features.

Proper orthogonal decomposition algorithm finds a set of optimal orthogonal bases${φi|i=1,2,…n}$ such that the projection error of the fluctuation part of sample set ${f′(θk)}k=1M$ in the space spanned by $Φ=[φ1,φ2,…,φn]$ is minimized, that is to say, $Q=1M∑k=1M(‖f′(θk)−Φαk‖)$ is minimized. Here, $αk$ is the vector composed of orthogonal polynomial expansion coefficients.

Solving its eigenvalues and eigenvectors, M non-negative eigenvalues, $λi(i=1,2,…M),λ1≥λ2≥⋯λM$⁠, and corresponding eigenvectors, $bi(i=1,2,…M)$⁠, are obtained. The optimal orthogonal base $φi$ can be expressed by $φi=∑k=1Mbkif′(θk)$⁠.

Through POD, the original high-dimensional response is reduced to an s-dimensional response, effectively alleviating the challenge of constructing a surrogate model. In subsequent artificial neural network modeling, the output is no longer the original response $f(θ)$⁠, but a low-dimensional response vector $α$ composed of basis mode coefficients. It is noteworthy that the formula (3) facilitates a bidirectional process, enabling not only dimensionality reduction of the known $f(θ)$ to obtain $α$ but also the utilization of $α$to reconstruct the original output response $f(θ)$⁠.

After obtaining reduced-dimension responses by POD, it is necessary to construct a surrogate model between the model uncertainty parameters and the reduced-dimension responses. The surrogate model that can be used here is not unique. The primary reason for selecting artificial neural networks is their ability to fit any measurable function with arbitrary accuracy and their excellent scalability with the size of the dataset [19].

The widely used error backpropagation algorithm is utilized to continuously adjust all the weights and biases of the entire network to minimize the loss function.

This surrogate model facilitates the prediction of base mode vectors $αnew$ for new model parameters θ_new, subsequently enabling the retrieval of the original output response $f(θnew)$ via POD.

To mitigate the dimensionality of the input uncertainty space, sensitivity analysis techniques are required to identify the factors exerting substantial influence on the output. Scholars have devised numerous global or local sensitivity analysis methodologies based on derivatives, variance, and other metrics. However, for high-dimensional outputs, it is imperative to analyze the collective impact of model parameters on the multidimensional correlated response, rather than conducting sensitivity analysis independently for each output dimension. To achieve this, this paper utilizes a global sensitivity analysis approach proposed by Gamboa et al. which decomposes the covariance of multi-output variables [21].

The main effect denotes the impact of the individual variation of an input variable on the output variance, whereas the total effect encompasses both the individual influence and the coupled effect of this input variable with other variables. As depicted in Eq. (13), the sensitivity index (comprising the main and total effects) of any parameter is obtained by weighting the sensitivity index of this parameter for each output dimension with the variance of that output component. For each output dimension, Monte Carlo sampling can be utilized to estimate the main and total effects of the input variables. In this study, a surrogate model is employed to substitute the original CFD simulation, facilitating efficient computation of large-scale samples.

To assess the efficacy of the proposed approach, this paper examines the influence of uncertainty in the coefficients of the SA turbulence model on the prediction of airfoil and wing aerodynamics. The SA model is extensively utilized in aerospace and related fields [22]. The computation presumes fully turbulent flow and disregards the transition term in the original model, resulting in nine constants denoted as: $cb1,σ,cb2,κ,cw2,cw3,cv1,ct3,ct4$⁠. It is assumed that the model coefficients are subject to epistemic uncertainty, which can be characterized using probabilistic methods from a credibility perspective. The model parameters are assumed to follow uniform distributions with intervals specified in Table 1, with parameter ranges referenced from pertinent literature [23]. It is worth noting that the mathematical representation of model uncertainty parameters, including distribution types and parameters, necessitates thorough consultation with model developers. The present study primarily focuses on the propagation of uncertainty given a mathematical description of the uncertain parameters, rather than on the precise quantification of uncertainty for the model parameters. Numerous investigations have been conducted on the quantification of uncertainty in SA turbulence model parameters, primarily employing polynomial chaos methods, with an emphasis on overall aerodynamic outputs such as overall forces [23,24].

The two numerical examples investigated are the distribution of wall friction coefficient in the low-speed flow over a NACA0012 airfoil and the distribution of wall pressure coefficient in the transonic flow over an M6 wing. These two cases will be presented separately as follows:

The NACA0012 airfoil, with a symmetric profile and 12% thickness, is considered under computational condition of $M∞=0.15,α=5.0∘,T∞=288.15K,Re∞=6×106$⁠. The computational grid employed is depicted in Fig. 4. The distribution of wall friction coefficient along streamwise coordinate obtained utilizing the standard parameters of the SA model is presented in Fig. 5, revealing relatively large changes at the leading edge. Accurately capturing these local, sharp changes in the friction coefficient pose considerable challenges for the prediction model. The output in this case comprises friction coefficients for every grid point on the entire airfoil surface, yielding a 64-dimensional result.

The transonic flow over the M6 wing serves as a benchmark test case for assessing transonic flow solvers. The wing features a root chord of approximately 0.8 m and a half-span of approximately 1.2 m, with a symmetric airfoil section. The computational condition is set at$M∞=0.8395,α=3.06∘,T∞=255.56K,Re∞=1.172×107$⁠. The computational grid employed in the calculation is displayed in Fig. 6. Under this condition, a $λ$-shaped shock structure develops on the upper surface of the wing, as shown in Fig. 7, posing substantial challenges for the prediction model. The output of this case constitutes a vector comprising pressure coefficients of all grid points covering the entire wing surface, with a dimensionality of 29,684.

The sample data in this paper are generated by in-house unstructured grid solver Flowstar [25], which is founded on a cell-centered finite volume methodology and is adept at handling diverse element types, including hexahedra, tetrahedra, prisms, pyramids, and other polyhedra generated via geometrical multigrid techniques. Second-order accuracy in space is attained through linear reconstruction within cells. The vertex-based Green-Gauss approach [26] is employed for gradient computations to uphold accuracy and robustness. To mitigate oscillations in regions of high gradients, Venkatakrishnan's limiter [27] is utilized. The Roe scheme is engaged for inviscid flux computations. A first-order backward Euler time-differencing scheme with local time stepping is implemented to approximate steady-state, facilitating convergence. The flux Jacobian is derived from a first-order upwind scheme, with the divided convective flux Jacobian composed of the convective flux Jacobian and its spectral radius. The viscous flux Jacobian is approximated using its spectral radius.

The prediction accuracy of the model can be assessed through the evaluation of the prediction error of the output vector $f$⁠. For testing sample, the error vector is defined as $ε=fmodel−fCFD$⁠, serving for evaluating the accuracy of the model.

where n is the dimension of the output vector, ref is the range of the sample response, $ref=max(fCFD)−min(fCFD)$

To objectively assess the prediction model's generalization capabilities, Latin hypercube sampling was conducted twice separately within the 9-dimensional input space, generating 110 training samples and 10 testing samples. The training samples facilitated the development of a prediction model capable of forecasting friction coefficient distributions on the airfoil surface for any model parameters. The testing samples were utilized to validate the prediction model's accuracy. Figure 8 illustrates the projection distribution of all samples within the two-dimensional parameter space of σ − κ.

Each sample input was processed by the CFD solver to derive the corresponding distribution of wall friction coefficients, constructing a comprehensive sample dataset. The snapshot matrix, composed of all training samples, undergoes POD decomposition. Basis functions are then truncated according to the 99.9% generalized energy criterion, leading to the retention of eight basis functions. Figure 9 illustrates the relationship between the energy proportion and the number of POD modes.

As per the methodology outlined previously, the outputs of all training samples were initially projected onto the basis function space, reducing the 64-dimensional friction coefficient vector to an eight-dimensional basis function coefficient vector. Subsequently, an artificial neural network was leveraged to develop a prediction model linking the nine-dimensional uncertain parameters and the eight-dimensional basis function coefficient vector. In this instance, the artificial neural network employed a 2-layer hidden layer architecture with 50 neurons in each layer. Once model training was completed, the basis function coefficient vector for test sample can be predicted via this artificial neural network. The corresponding distribution of wall friction coefficients was then derived using the bidirectional representation of POD. By comparing the model's predicted outcomes with CFD calculation results, the prediction error of the model can be assessed. The prediction error was computed for all testing samples, and the resulting error distribution is depicted in Fig. 10. Across all testing samples, it can be noted that $ε1$ remains below 0.035%, $ε2$ remains below 0.06%, and $εinf$ remains below 0.2%, maintaining at a fairly low level.

Two test samples exhibiting relatively high prediction errors, specifically samples No. 1 and 3, were chosen for detailed comparison. As shown in Fig. 11, the predictions from the model demonstrated very close agreement with the CFD calculations, even at the leading edge, where sharp changes occur. Consequently, this prediction model is suitable for conducting large-scale Monte Carlo sampling to investigate the uncertainty of the friction coefficient distribution. The time consumed by the model for prediction is negligible in comparison to CFD simulations.

Finally, the uncertainty of friction coefficient is analyzed. For practical engineering problems, there is no true value for uncertainty. In this paper, the random sampling results based on CFD simulations are used as the reference. Considering that the random sampling results are significantly related to the sample size, the statistical information under different random sample sizes is analyzed first, as shown in Fig. 12. All the random samples in the figure are obtained through CFD simulations, not by prediction models. When increasing the number of samples from 110 to 1100, the mean value of the friction coefficient remains relatively stable, but the standard deviation exhibits notable variations. Notably, the standard deviation obtained using 550 samples is comparable to that obtained using 1100 samples. Therefore, for the purpose of comparison, we consider the statistical results obtained using 1100 samples as the reference.

A large-scale random sampling analysis of the uncertainty of the friction coefficient was conducted with the prediction model. Due to the negligible prediction cost of the prediction model, Latin hypercube sampling was used to obtain 10⁶ samples. The mean value and standard deviation of the friction coefficient are shown in Fig. 13. The prediction results of this model are in good agreement with the direct statistical results of 1100 CFD samples. Notably, the uncertainty of the friction coefficient is relatively significant across the entire airfoil, particularly at the front part of the upper surface.

The relative significance of the model parameters on the distribution of friction coefficients can be ascertained through the sensitivity analysis approach based on covariance decomposition, as presented in the Table 2. κ emerges as the most influential parameter, which can be expected from a physical standpoint. In the constant Reynolds stress layer of plane shear turbulence, the eddy viscosity coefficient is directly proportional to κ, i.e., ν_t = κu_τd. The eddy viscosity coefficient characterizes the relationship between Reynolds stress and mean flow, serving a role analogous to molecular viscosity in RANS simulations. Hence, κ is intimately linked with the viscosity of the flow, and the friction coefficient is directly correlated with the flow viscosity.

Similar to the first case, Latin hypercube sampling was conducted twice separately within the 9-dimensional input space to generate 110 training samples and 10 testing samples. The training samples facilitated the development of a prediction model capable of forecasting pressure distributions on the wing surface for any given SA model parameters, while the testing samples were utilized to validate the prediction accuracy. Each sample input was processed by the CFD solver to derive the corresponding wall pressure distribution, constructing a comprehensive sample dataset. In principle, an artificial neural network or multi-output Gaussian model could be directly employed to establish a prediction model between the 9-dimensional uncertain parameters and the 29,684-dimensional output. However, due to the significant number of model parameters to be estimated and the prolonged training time, a flow field reduction technique is necessary to reduce the output dimensionality.

All training sample snapshots undergo POD decomposition, and the basis functions are truncated according to the 99.9% generalized energy criterion, yielding nine basis functions. The outputs of all training samples are then projected onto the basis function space, reducing the dimensionality of output from 29,684 to 9. An artificial neural network is leveraged to develop a prediction model that maps the 9-dimensional uncertain parameters to the nine-dimensional basis function coefficient vector. Specifically, the artificial neural network utilizes a 2-layer hidden layer architecture with 100 neurons in each layer. After completing model training, the basis function coefficient vector for any test sample can be predicted via the artificial neural network. Subsequently, the corresponding wall pressure vector is derived through the bidirectional representation of POD. By comparing the predicted outcomes of the entire model with CFD calculation results, the prediction error can be quantified. The resulting error distribution is depicted in Fig. 14. Across all testing samples, it can be noted that $ε1$remains below 0.0015%, $ε2$ remains below 0.008%, and $εinf$ remains below 0.4%, maintaining at a fairly low level.

To conduct a more intuitive validation of the prediction accuracy, the fourth test sample was chosen for detailed comparison. Pressure coefficient contours predicted by the model and through CFD simulation, along with comprehensive comparisons of six cross section spanning from the root to the tip of the wing (cross section locations illustrated in Fig. 7) are provided in Figs. 15 and 16. In Fig. 16, the sample mean was subtracted to enhance visibility of the discrepancies between the two. The displayed results indicate a high degree of concordance between the prediction model and CFD calculations, even at the complicated shock wave area. Consequently, this prediction model is suitable for large-scale Monte Carlo sampling to investigate uncertainty information related to the pressure coefficient distribution. The time consumed by prediction model is insignificant in comparison to CFD simulations.

Ultimately, Latin hypercube sampling was employed again within the nine-dimensional input space to generate 10⁶ samples. These inputs were sequentially fed into the fast prediction model to acquire the corresponding pressure coefficient vectors. The contours of mean value and standard deviation of the pressure coefficient are shown in Fig. 17. Notably, the uncertainty of the pressure coefficient is concentrated at the $λ$-shaped shock on the upper surface, whereas the uncertainty at other locations is minimal. This observation aligns with expectations, given that the prediction of shocks is highly sensitive to turbulence model parameters. Similarly, we use random sampling results based on CFD simulations as a reference. Due to the significantly higher computational cost of this case compared to the previous one, we did not perform a rigorous analysis of the impact of sample size. However, a comparison of standard deviation distributions across the six sections in Fig. 18 reveals that the method proposed in this paper is highly similar to the results obtained from 550 CFD random samples. Given that the prediction model was built using only 110 CFD simulation samples, the prediction cost of the prediction model is significantly lower compared to random sampling methods.

Through the application of the multi-output global sensitivity analysis approach grounded in covariance decomposition, the relative significance of model parameters on the pressure coefficient distribution can be ascertained, as presented in the Table 3. Again, κ emerges as the most influential parameter. Following the same analysis as the first case, κ is intimately linked with the viscosity of the flow, which significantly impacts the kinetic energy loss of fluid microgroups within the boundary layer and the capacity to resist adverse pressure gradients, ultimately influencing the prediction of shock wave positions.

This paper addresses the need for quantifying uncertainties in multidimensional correlated responses within fluid flow fields. By employing a flow field dimension reduction approach based on proper orthogonal decomposition and a modeling methodology using artificial neural networks, we achieve rapid prediction of flow field variables under uncertain model parameters. Utilizing this fast prediction model, Monte Carlo method was utilized to analyze the uncertainty of flow field variables. Furthermore, a multi-output global sensitivity analysis approach based on covariance decomposition is adopted to identify significant model parameters that impact the multidimensional correlated responses of the flow field. The results indicate that

1. Proper orthogonal decomposition enables the identification of key basis functions representing flow structures using a limited number of sample snapshots. This effectively reduces the dimensionality of flow field correlated responses from hundreds or thousands of dimensions to an about ten dimensional basis function coefficient vector, significantly reducing the output dimensionality of uncertainty propagation analysis. Additionally, the bidirectional representation of flow field variables and basis function coefficients through proper orthogonal decomposition, combined with artificial neural network modeling, enables rapid prediction of flow field variables under any model parameters. The two examples used in this study demonstrate that the dimensionless maximum error of fast model prediction is within 0.4%, with accurate predictions of flow field variables in rapidly varying regions, providing a reliable prediction model for uncertainty propagation analysis.

2. The multi-output global sensitivity analysis approach based on covariance decomposition identifies the parameters that have a significant overall impact on flow field variables, overcoming the limitations of individual output dimension analysis methods. For the two examples in this study, κ is found to be the most critical factor, consistent with physical intuition.

The sample data utilized in this study were acquired at a consistent level of accuracy. Currently, multifidelity modeling approaches have been widely applied in industrial optimization and design fields. The introduction of low-fidelity sample data demonstrates lower overall costs compared to solely using high-fidelity sample data. Therefore, future efforts will focus on developing rapid and accurate prediction models for fluid flow fields, leveraging multifidelity sample data and adaptive parameter space sampling techniques to enhance the efficiency of the prediction models and improve the effectiveness of uncertainty propagation analysis.

• NSAF (Grant No. U2230208; Funder ID: 10.13039/501100001809).

• National Numerical Wind Tunnel Project (Funder ID: 10.13039/501100018618).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.