Abstract
Airfoil geometric uncertainty can generate aerodynamic characteristics fluctuations. Uncertainty quantification is applied to compute its impact on the aerodynamic characteristics. In addition, the contribution of each uncertainty variable to aerodynamic characteristics should be computed by the uncertainty sensitivity analysis. In the paper, Sobol’s analysis is used for uncertainty sensitivity analysis and a nonintrusive polynomial chaos method is used for uncertainty quantification and Sobol’s analysis. It is difficult to describe geometric uncertainty because it needs a lot of input parameters. In order to alleviate the contradiction between the variable dimension and computational cost, a principal component analysis is introduced to describe geometric uncertainty of airfoil. Through this technique, the number of input uncertainty variables can be reduced and typical global deformation modes can be obtained. By uncertainty quantification, we can learn that the flow characteristics of shock wave and boundary layer separation are sensitive to the geometric uncertainty in transonic region, which is the main reason that transonic drag is sensitive to the geometric uncertainty. The sensitivity analysis shows that the model can be simplified by eliminating unimportant geometric modes. Moreover, which are the most important geometric modes to transonic aerodynamics can be learnt. This is very helpful for airfoil design.
1. Introduction
There inherently exist a vast of uncertainties in the practical aircraft design and application, thus causing fluctuations of aircraft performance. Therefore, it is important to take these uncertainties into account at the beginning of aircraft design. With the development of computer technology, computational fluid dynamics (CFD) technology has been widely used to solve problems in aerodynamic mechanics. The traditional CFD simulation is deterministic. However, a variety of uncertainties are inevitable in CFD simulation with the increasing complexity of the fluid problem in reality, leading to the mismatch between the CFD simulation results and the actual results [1, 2].
The sources and classifications of uncertainty in CFD were described in [3, 4]. Several uncertainty quantification (UQ) methods have been used in CFD, including Monte Carlo simulation (MCS), sensitivity analysis (SA), moment methods, and polynomial chaos expansions (PCE) in [5]. MCS is a statistical method, which needs a large number of samples to accurately analyze the uncertainty quantification. SA and moment methods are suitable to solve the problem of small parameter uncertainty or linear model. Recently, PCE has been widely applied to UQ of fluid problems. PCE methods can be divided into intrusive and nonintrusive methods according to the coupling ways with CFD solvers. In general, an intrusive approach calculates the unknown polynomial coefficients by projecting the resulting equations into basis functions for different modes, and it requires modifying the CFD codes, which may be difficult and time-consuming for complex problems such as N-S simulation. Mathelin et al. [6] used the intrusive PCE to research the quasi one-dimensional duct flow uncertainty propagation problems. Xiu and Karniadakis [7] extended the original polynomial methods to the general polynomial chaos (gPC) and applied it to the incompressible channel flow and the flow around a cylinder. To overcome the shortcomings of intrusive polynomial chaos, nonintrusive polynomial chaos (NIPC) has been developed. CFD is regarded as a black box without changing the CFD program code in nonintrusive methods. There are two different methods to build a PCE model: projection method and regression method. The projection method is based on the orthogonality of PCE to evaluate the coefficients by a multidimensional integral. There are two approaches to evaluate this integral: random sampling method and Gaussian quadrature-based method. Random sampling method uses MCS to evaluate the unknown coefficients, but its convergence rate is low. Deterministic sampling method uses the quadrature to evaluate the unknown coefficients. The quadrature-based methods are more efficient than random sampling methods in low-dimensional problems. However, they are inefficient for relatively high-dimensional problems because of the exponential rising in quadrature points with the increasing dimensions. The other nonintrusive collocation method evaluates PCE using regression. Hosder et al. [8] developed this regression method to NIPC for stochastic computational fluid dynamics. UQ based on PCE and its applications in fluid mechanics were comprehensively reviewed in [9, 10].
In addition to UQ, global sensitivity analysis (GSA) plays an important role in quantifying the relative importance of each uncertainty source. By means of this technique, uncertainties can be systemically studied to measure their effects on the system outputs so as to screen out the uncertainties with negligible contributions to simplify the problem. A comprehensive review of uncertainty sensitivity analysis can be found in [11]. The variance-based method, also called Sobol’s analysis [12, 13], is one of the most popular practices in many disciplines. Sobol’s analysis can measure the relative importance of one input variable to the model’s output by computing the partial variance. There are many methods to compute Sobol’s indices, such as Fourier Amplitude Sensitivity Test (FAST) and sampling-based method [14, 15]. Sudret [16] introduced PCE to build surrogate models that allow one to compute Sobol’s indices analytically as postprocessing of PCE coefficients.
In stochastic aerodynamic analysis, most research on UQ studies the flight condition parameter uncertainty based on NIPC methods. Loeven et al. [17] conducted a subsonic aerodynamic analysis around a NACA0012 airfoil with an uncertain free stream velocity using a commercial flow solver. From their study, an uncertain free stream velocity leads to the highest variation in pressure on the upper surface near the leading edge. Simon et al. [18] focused on the transonic stochastic response of two-dimensional airfoil to parameter uncertainty (Ma and ) using gPC. Two kinds of nonlinearities are critical to transonic aerodynamics in their study: the leeward shock movement characteristics and boundary layer separation on the aft part of airfoil downstream the shock. Chassaing and Lucor [19] conducted a stochastic investigation of flows about NACA0012 airfoil at transonic speeds. Liu et al. [20] conducted a stochastic fluid analysis on a 3D wind blade case considering the wind speed as an uncertain parameter. From their study, when the flow separation appears, the separation vortex region corresponds to the maximum variation area, and the maximum variation extends to the trailing edge even to the whole suction side. Moreover, Sobol’s analysis was also introduced to stochastic aerodynamic analysis. Simon et al. [18] applied Sobol’s analysis to the stochastic transonic aerodynamic analysis considering the uncertainty of Mach number and angle of attack. Hosder and Bettis [21] conducted an uncertainty and sensitivity analysis for reentry flows with inherent and model-form uncertainties. However, rare research can be found on stochastic aerodynamics analysis considering geometric uncertainty. The geometric uncertainty on aerodynamic surfaces resulting from manufacturing errors has significant effect on the aerodynamic performance. Due to the high cost of precise surface manufacturing techniques, it is often impractical to remove the impact of these geometric variations by improving the manufacturing tolerance. In other words, the geometric uncertainty resulting from manufacturing errors is unavoidable. Therefore, it is important to conduct an UQ for aerodynamics evaluation considering this geometric uncertainty. Moreover, a GSA for aerodynamics evaluation considering geometric uncertainty is also necessary, because it can obtain important geometric variation parameter to aerodynamics by sensitivity analysis. It is difficult to describe geometric variation because it needs a lot of input parameters. The high dimensions of input uncertainty can lead to the high computational cost to build PCE model, which makes it more difficult to conduct a stochastic aerodynamic analysis considering geometric uncertainty.
In this paper, we conduct UQ and GSA of transonic aerodynamics of airfoil RAE2822 considering geometric uncertainty. The principle component analysis (PCA) technique combining with an airfoil parameterization method is used to describe the airfoil shape deformation and obtain the typical deformation modes. Then, a stochastic PCE model is built for UQ and GSA. The paper is structured as follows. Section 2 introduces UQ and the variance-based GSA (Sobol’s analysis). In Section 3, two methods, MCS and PCE, are introduced to UQ and GSA in detail. In Section 4, UQ and GSA of aerodynamics are conducted. Section 5 outlines several useful conclusions.
2. Uncertainty Quantification and Global Sensitivity Analysis
The aim of UQ is to compute the influence of input parameters on the model’s output quantitatively. The input uncertainty may be expressed in a number of ways, for example, interval bounds or probability density functions. In the paper, the uncertainty is quantified by the probabilistic analysis and the robustness (mean and variance) of outputs is focused on. The mean and variance are shown by the following equation:
Although the influence of the input uncertainty on output variables can be identified by UQ, it cannot obtain the contribution of each uncertainty. In this case, A GSA is applied to study how the output of a model can be apportioned to different sources of input uncertainty. Since most of the sensitivity analyses are based on the partial derivate, this approach is sometimes called local sensitivity. However, the fatal limitation of a derivative-based approach is that when the model is uncertain is unwarranted. In other words, derivatives are only informative at the base point where they are computed and do not provide an exploration of the rest of the uncertainty space. The alternative uncertainty sensitivity analysis is GSA. The main goal of GSA is to demonstrate the relative importance of each input uncertainty among the overall uncertainty in the output. It is often useful to rank the contribution of each uncertainty to the overall uncertainty in the output. A detailed review of a variety of GSA methods can be found in [15]. In the current study, a variance-based global sensitivity method is used to accomplish this task. The method is based on Sobol’s variance decomposition, by which -function can be uniquely decomposed into functional terms of increasing dimensions:where , , and . Adding variance to both sides of (2), we can obtain the following equation:where . The first-order partial variance can be regarded as the average reduction of model output variance resulting from fixed ; that is, measures the individual contribution of to the total variance . The larger is, the more reduction of output variance can be obtained by reducing the uncertainty of . The second-order partial variance quantifies the interaction effect between and . Similar interpretations can be given to the higher-order partial variances. Dividing both sides of (3) by , we obtain
are Sobol’s indices, which measure the relative contribution to the total variance. are the main sensitivity indices. And measure the interaction effects.
Another commonly used measurement is the total partial variance , which is defined as the summation of all terms in (6) with subscripts including . In other words, incorporates both the individual effect of and its interaction effects with all the other input variables.
Usually is used to select important variables, while is more suitable to screen noninfluential variables.
3. The Nonintrusive Polynomial Chaos
There are many uncertainty analysis methods available for UQ [22, 23]. In the paper, The NIPC method is used for the current study. The PCE is a stochastic method based on the spectral representation of uncertainty. According to the spectral representation, the random function can be decomposed into deterministic and stochastic components. For example, a random variable () can be represented by the following equation:where is the deterministic component and is the random basis function corresponding to the th mode. From (8), the random variable is the function of deterministic independent variable vector and the -dimensional standard random variable vector . The polynomial chaos expansion given by (8) contains an infinite number of terms. In a practical computational context, the terms of PCE can be truncated by both order and dimension . The number of terms is finite, which is given by the following equation:where is the number of random dimensions and is the order of polynomial chaos. When the input uncertainty obeys Gauss distribution, the basis function is the multidimensional Hermite polynomial.
There are two ways to solve the coefficients of (8): intrusive and nonintrusive methods. The intrusive method computes the unknown polynomial coefficients by projecting the resulting equations into basis functions for different modes. It requires the modification of the deterministic code, which may be difficult, expensive, and time-consuming for complex computational problems such as complex N-S equations. Alternatively, the nonintrusive method treats the CFD as a black box without changing the program code when propagating uncertainty. Now, we pay close attention to how to solve the coefficients in nonintrusive methods. There are two methods: projection method and regression method.
The main idea of projection method is to use Galerkin projection to solve coefficients. Equation (8) can be transformed to (10) by inner product: represents the inner product, which can be expressed by the following formula:For Gaussian random variable, the basis functions are Hermite orthogonal polynomials. Because of orthogonality, (10) can be transformed to the following equation:And then the following can be derived:
The sampling approaches can be divided into random sampling method and deterministic sampling method. Random sampling method uses MCS, LHS, and other methods to compute the projection integrals. However, its convergence rate is low. Deterministic sampling method uses the quadrature for the numerical evaluation of the unknown coefficients. Using -dimensional Gauss-Hermite quadrature, with points in each dimension, we can compute the unknown coefficients by the following equation:
In (14), one-dimensional integral is expanded into a high-dimensional form by tensor product, the calculation times of which require points for th-order polynomial chaos. For low-dimensional problems, the sampling efficiency of deterministic sampling method has been greatly improved compared with random sampling method. However, the calculation times grow exponentially with the increasing dimensions.
The regression method is another approach to compute the coefficients using a linear equation system based on a selected set of points. A linear system of equations can be obtained bywhere is given by . The accuracy of regression method depends on the selection of sample points in design space. According to [16], the optimal design is given by the roots of Hermite polynomial, in which the optimal sample points are . The sample point of NIPC with projection method is , which grows exponentially with the increasing number of input dimensions. Compared with projection based method, the sample points can be reduced by regression based method. The sample points are , which gives a better approximation at each polynomial degree in [8]. To study how the number of samples influences the analysis results, the oversampling ratio is defined as
After computing the unknown coefficients, the approximate PCE model is built, and then the model is used to compute the first two statistical moments analytically as follows:
In order to compute Sobol’s indices, Sobol’s decomposition based on PCE should be derived. The detailed derivation process is in [16]. The important equations are shown:where is the th Hermite polynomial. Let us define the set , where the indices are nonzero:
Sobol’s analysis of PCE is shown by (19). Combining (4)–(7), Sobol’s indices are obtained.
4. Uncertainty Quantification and Sensitivity Analysis for Transonic Aerodynamics
4.1. The Description of Geometric Uncertainties
In engineering, the geometric uncertainty on aerodynamic surfaces resulting from manufacturing errors or wearing is unavoidable, which has significant influence on the aerodynamic performance. Therefore, UQ and GSA considering geometric uncertainty to aerodynamics are conducted in this section. To carry out this work, the first step is to describe the geometric uncertainty in the computing environment.
In [24, 25], the main geometric variation modes are obtained by principle component analysis (PCA) with a large amount of geometric statistical manufacturing error data of airfoil. The main variation modes based on the statistical manufacturing error data are used to achieve the geometric variation. The description of the geometric uncertainty of airfoil is shown in the following equation: where is the nominal geometry; is the average geometric variation; is the geometric mode shape; is the number of mode shapes used to represent the variation in geometry. The geometric mode can be computed by PCA based on manufacturing samples. is the th singular value of the measurement snapshot matrix, which represents the geometric variability attributable to the th mode. is a random parameter which obeys the standard normal distribution; thus, the product is the stochastic contribution of th mode.
It is difficult to describe this geometric variation in the computing environment. In [26], a Gaussian random process simulation is used to obtain the geometric data. Then PCA is used to obtain main geometric modes. In addition, a variety of parametric methods have been employed to describe the geometric variation in aerodynamic design so far, like PARSEC-11 geometry parameterization, class-shape function transformation (CST) method, free-form deformation (FFD) method, Hicks-Henne bump functions, and so forth [27–31]. By changing the parameter of these parametric methods, the geometric variation is realized. Generally, these parametrization methods need a lot of parameters to represent the airfoil shape. To reduce the dimensions of the variables, the PCA technology combined with airfoil parameterization is used to describe the geometric variation used in the paper. Firstly, a parametric method is used to generate a set of sample data. Then the PCA is used to obtain the main deformation modes based on the generated sample data. In this paper, we use a parameterized representation of the airfoil RAE2822 by CST method with 24 parameters. The data of the measurement points on the airfoil surface are obtained by random perturbation of CST parameters. In this way, the PCA based sample data is conducted. Figure 1 shows the eigenvalue with the number of modes, which represents the geometric variation attributable to each mode. The eigenvalue of one mode is smaller; the proportion of the mode to geometric variation is smaller. It can be observed that the proportion of the first 12 modes is 99.32%, which means that they can be used to express the geometric variation well. Through the PCA, the number of input parameters is reduced from 24 to 12. The first 12 modes obtained by PCA are shown in Figure 2. The 12 modes include 6 modes of the upper surface (mode 2, mode 3, mode 5, mode 7, mode 9, and mode 10) and 6 modes of the lower surface (mode 1, mode 4, mode 6, mode 8, mode 11, and mode 12). It can be seen that the lower-order modes in these modes are global geometric deformation modes and present some typical geometric deformation. Specifically, mode 1 and mode 2 are the scale modes in the thickness direction; mode 3 and mode 4 are translation modes of the maximum thickness in the axial direction; and mode 5 and mode 6 are the extrusion modes of the upper surface. So, the number of parameters used to describe the geometric variation is reduced to a large extent and the typical deformation modes are obtained through the PCA technique.
4.2. Deterministic Flow Computations (CFD Simulations)
The deterministic steady-state flow solutions are computed by the Reynolds-Averaged Navier-Stokes (RANS) equations. The two-dimensional Navier-Stokes equations can be written in the integral form aswhere is the control volume; is the boundary of control volume; is the outer normal vector of the control volume boundary; is the vector of conservative variables; is the inviscid flux; and is the viscous flux.
For the cell-centered finite volume method, the computational domain is divided into nonoverlapping control volumes that completely cover the domain. The interface variables are derived from the average values of the grid cells in order to calculate the flux of control volumes. Through spatial discretization, the equations of the integral form are transformed to the linear ordinary differential equations, and the flow variables can be obtained by the time marching method. We adopt AUSM scheme to evaluate the numerical flux which can also be evaluated by the upwind scheme [32]. The geometry of airfoil is changing in the process of UQ, and radial basis functions (RBF) are used for mesh deformation [33]. Figure 3 shows the computing grids of RAE2822 and the reliability verification of the CFD program.
(a) Global grid
(b) Local grid
(c) The comparison of
4.3. Uncertainty Quantification of Aerodynamics
UQ of aerodynamics is conducted in this section. The uncertain variables are in (20), which obey Gaussian distribution and are independent from each other. The flow state is selected in a transonic region (Ma = 0.73, °, and ); NIPC and MCS are used to quantify uncertainty. MCS is used to aim at validating the accuracy and efficiency of NIPC. The computational cost of MCS is selected as 5000 samples. Figure 4 shows the range of geometric uncertainty in the paper. It can be seen that the range of fluctuation reaches 0.006 near the maximum thickness, which is 10% of the maximum thickness. In this paper, the fluctuation range of the geometry should be greater than the fluctuation caused by the manufacturing error.
4.3.1. The Difficulty of NIPC for the High-Dimensional Problem
For a high-dimensional problem, the number of PCE terms is very high. Now, we study the number of terms with increasing for the 12-dimensional problem. From Table 1, it can be seen that the number of terms increased rapidly with increasing for the 12-dimensional problem.
It is mentioned in (14) that the projection based NIPC method needs samples for UQ, which grows exponentially with the increasing number of input dimensions. Table 2 shows the number of samples of projection based NIPC method with for the 12-dimensional problem. It can be seen that the computational cost is too high and unacceptable.
The regression based NIPC method adopts the linear regression method to calculate the coefficients of PCE. To study how the number of samples influences the analysis results, the oversampling ratio is defined as (16).
Compared with the projection based NIPC method, the computational cost of the regression based NIPC method is much lower. Therefore, the regression based NIPC method is used for the stochastic aerodynamic analysis considering geometric uncertainty.
4.3.2. The Verification of NIPC for the UQ considering Geometric Uncertainty
The verification of the regression based NIPC method for the UQ considering geometric uncertainty is conducted.
Firstly, we study the convergence of NIPC by increasing the polynomial chaos order . Figures 5–7 show the standard deviation of aerodynamics with the order of NIPC and the comparison with that of MCS. It can be seen that when or 3, the results of NIPC can coincide with those of MCS. However, there is a little difference between NIPC and MCS. It can also be observed that the results of NIPC are wrong when . We think that the number of terms of PCE is 1820 for and , which makes it difficult to compute the coefficients correctly by linear regression method. This leads to the wrong results for high-order accuracy.
Moreover, the convergence of NIPC by increasing the oversampling ratio is also studied. Figures 8–10 show the standard deviation of aerodynamics with of NIPC and the comparison with that of MCS. When , the results of NIPC can coincide well with those of MCS. Through comparison, it can be observed that the results of NIPC with and are even better than those of NIPC with and .
Lastly, the computational costs are compared with MCS. Table 3 shows the number of samples to UQ by NIPC and MCS. It can be seen that the number of samples is 364 with and , which is much lower than that of MCS.
4.3.3. The Results of UQ for Transonic Aerodynamics
Figure 11 shows the error bars of and distribution by UQ. It can be observed that the aerodynamic loads are sensitive to the geometric variation of the upper surface. In transonic region, the location of shock fluctuates with the change of the geometric shapes, which can result in intensive fluctuations of the flow variables. We call the region in which the standard variations become much larger suddenly shock disturbance region. It can also be seen that the skin-friction behavior displays discrepancies from the distribution. The variation of is similar to that of before and in the shock disturbance region. However, the difference occurs in the downstream region of the shock disturbance region; the standard deviation keeps larger magnitudes below the shock disturbance region and then gradually decreases. This indicates that not only the shock wave but also the boundary layer separation is sensitive to geometric uncertainty in the transonic region.
4.4. Global Sensitivity Analysis of Aerodynamics
The important deformation modes can be found out through GSA, which can provide a useful reference for subsequent design work and manufacturing process. In the paper, the NIPC method (; ) is used for Sobol’s analysis. The aerodynamic load distribution and aerodynamic coefficients are included for Sobol’s analysis. Figures 12 and 13 show the results of Sobol’s analysis of the pressure and friction coefficient distribution. The partial standard deviation, total standard deviation, and standard deviation caused by the coupling effect for each mode are given in Figures 12 and 13. It can be seen that the deformation modes of the upper surface are much more important than those of lower surface for transonic aerodynamic loads ( and ). Specifically, the first 3 modes are important for aerodynamic load distribution. For distribution, the main variation of modes 2, 3, and 5 is in shock disturbance region. For distribution, the main variation of modes 2 and 5 is in both shock disturbance region and boundary layer separation region, while the main variation of mode 3 is only in shock disturbance region. In addition, it can also be observed that the coupling effect among geometric deformation modes is weak in Figures 12 and 13. In order to visually observe the effects of modes, the pressure coefficient and friction coefficient of error bars are given in Figures 14 and 15. From Section 4.1, we can observe that mode 1 and mode 2 are the scale modes in the thickness direction; mode 3 and mode 4 are translation modes of the maximum thickness in the axial direction; and mode 5 is the extrusion mode of the upper surface. We can draw the following conclusions: the scale mode and extrusion mode are important for the flow characteristic of shock wave and boundary layer separation, but the translation mode is only important for flow characteristic of shock wave.
The results of GSA for aerodynamic coefficients including lift coefficient (), drag coefficient , moment coefficient (), and lift-to-drag ratio () are shown in Figure 16. For , mode 1 and mode 2 are important, which indicate that the scale modes are important for . For and , modes 2, 3, and 5 are important; the first 3 deformation modes of the upper surface have bigger contribution. For , the scale mode of upper surface is most important. In addition, it can also be seen that the coupling effect among deformation modes is very weak.
(a) Sobol’s indices for
(b) Sobol’s indices for
(c) Sobol’s indices for
(d) Sobol’s indices for
From the results of GSA, it can be observed that modes 1, 2, 3, and 5 are important for aerodynamics. Therefore, the UQ of aerodynamics is conducted considering the 4 modes. Moreover, because the coupling effects among the modes are small from GSA, we adopt regression method with parameters and to build the PCE model without the interaction terms for UQ. The comparison of the UQ results is shown in Table 4. It can be observed that the UQ results are relatively small different from each other; and the number of sample points is reduced from 364 to 32. That is to say, the unimportant modes are eliminated and PCE model can be modified by GSA, which can reduce the computational cost in UQ and robust design.
5. Conclusion
In the paper, the regression based NIPC method was used for UQ and GSA for transonic aerodynamics with geometric uncertainty. It is difficult to describe the variation of geometric shapes, which always requires a large number of parameters. The high dimensions of input parameters can result in the high computational cost of NIPC. To alleviate the contradiction between the number of input dimensions and the computational cost, the PCA technique combined with CST method is used to describe the geometric uncertainty, which reduces the number of input variables. Moreover, typical geometric deformation modes can be obtained by the method.
From the results of UQ, the fluctuation of aerodynamic characteristics is mainly centered in the shock disturbance region and boundary layer interference region in transonic region. The characteristics of shock wave are sensitive to geometric uncertainty, and since the transonic drag mainly includes wave drag, the transonic drag is sensitive to geometric uncertainty. From the results of GSA of aerodynamic loads, it is shown that the first three deformation modes of the upper surface have the greatest impact on the aerodynamic characteristics. Specifically, the scale mode and extrusion mode are important for the flow characteristics of shock wave and boundary layer separation, but the translation mode is only important for the flow characteristics of boundary layer separation. Moreover, from the results of GSA of aerodynamic coefficients, we can observe that the scale modes (modes 1 and 2) are important for lift characteristics, while the contribution of the transition modes (modes 3 and 4) is negligible. The first 3 modes of the upper airfoil surface are important to drag and lift-to-drag characteristics. By GSA, we can eliminate unimportant geometric modes and simplify model for UQ and learn which typical modes are important for transonic aerodynamics, which is helpful for transonic airfoil design.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
The National Natural Science Foundation of China (no. 11572252), “111” Project of China (no. B17037), the National Science Fund for Excellent Young Scholars (no. 11622220), and the National Natural Youth Science Foundation of China (no. 11502209) supported this work.