A Two-Phase Monte Carlo Simulation/Non-Intrusive Polynomial Chaos (MSC/NIPC) Method for Quantification of Margins and Mixed Uncertainties (QMMU) in Flutter Analysis

A two-phase Monte Carlo Simulation/Non-intrusive Polynomial Chaos (MCS/NIPC) method for quantification of margins and mixed uncertainties (aleatory and epistemic uncertainties) is proposed in this paper for the flutter speed boundary analysis. Compared with the traditional MCS/MCS method which needs lots of numerical simulations, the MCS/NIPC method can reduce the computational cost without losing accuracy due to the use of point collocation non-intrusive polynomial chaos in the inner loop. Based on the results of uncertainty quantification, a novel practical quantification of margins and mixed uncertainties (QMMU) framework is proposed considering both aleatory and epistemic uncertainties such as the parametric uncertainties in complicated models. A two-dimensional airfoil and a three-dimensional benchmark wing, the AGARD 445.6 wing, are both employed to illustrate the practical application of the proposed methods for flutter speed computation in the presence of mixed uncertainties arising from the material properties and flight conditions. Within the analysis, aleatory and mixed uncertainty quantification are conducted to prove the effectiveness and efficiency of the MCS/NIPC method. Then the QMMU metrics are accomplished for the flutter speed boundary analysis of the two airfoil models. The results demonstrate the potential of the proposed QMMU framework to analyze the stability of engineering systems.


I. INTRODUCTION
Flutter is the dynamic instability where the structure extracts kinetic energy from air and this energy cannot be dissipated by structural damping. It is an important branch in the field of aeroelasticity, which has received a significant amount of research since 1960s [1], [2]. An accurate estimation of flutter characteristics including flutter speed has been considered as a crucial part in the aeroelastic analysis of aircraft design. Flutter analyses of three-dimensional panels [3], [4], The associate editor coordinating the review of this manuscript and approving it for publication was Yu Liu . composite laminated panels [5], swept high-aspect-ratio wing models [6] and a typical five degrees of freedom wing [4] have been conducted both analytically and experimentally. In addition, the effects of the boundary conditions, material properties, temperature fields, and yaw flow angles on the flutter characteristics were investigated [7].
However, the aforementioned studies of flutter speed were performed under the assumptions of deterministic situations where all the structural and aerodynamic parameters were known. As a matter of fact, in real aeroelastic systems there are multiple sources of uncertainties including but not limited to the modeling-induced uncertainties [8], [9], numerical VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ uncertainties, and parametric uncertainties [10]. The above uncertainties have a significant impact on the flutter speed boundary calculation. With the development tendency of precise and meticulous design for aeroelastic system, there is a growing concern in understanding how the uncertainties that inevitably exist in numerical models affect the predicted flutter speed boundary. Pettit [11], Dai and Yang [12], and Beran et al. [13] reviewed the uncertainty quantification in aeroelasticity, including both probabilistic and non-probabilistic methods. Typically, uncertainties in the aeroelastic systems could be divided into two commonly recognized categories: aleatory uncertainty and epistemic uncertainty [14]. For aleatory uncertainty which represented by a probability distribution, the probabilistic approaches based on the Monte Carlo Simulation (MCS) [15], the polynomial chaos expansion [16], the stochastic perturbation [17], [18], and the stochastic operator method [19] are powerful tools in aeroelastic analysis of the model with uncertainties. Jia et al. [20] developed an uncertainty propagation analysis method which was based on an extended sparse grid technique and maximum entropy principle to fit accuracy of the probability density function of the system response. The Non-intrusive Polynomial Chaos (NIPC) derived from the polynomial chaos expansion which is proposed by Ghanem and Spanos [21], has been widely used in the uncertainty propagation in the aerospace area [8]- [10]. For epistemic uncertainty which stems from the insufficient information about the input parameters, the nonprobabilistic methods including the interval theory [22], the fuzzy theory [23], the evidence theory [24], and the structured singular value (µ) method [12], [23], have been prevalently employed to quantify the propagation of uncertainties in the aeroelastic analysis.
Till now, most of the work regarding uncertainties in the aeroelastic analysis concentrates on flutter characteristics in the presence of a single uncertainty, either aleatory uncertainty or epistemic uncertainty. It is desirable to consider the contributions of both aleatory and epistemic uncertainties simultaneously and to quantify their influences on the flutter speed. To represent the impacts of two types of uncertainty, the boundary of cumulative distribution function (CDF) of the system response quantities called Probability box (P-box) [25] was utilized. To efficiently compute the probability bounds of system response function, new uncertainty propagation methods for problems with parameterized P-box were proposed [26], [27]. Meanwhile, the nested sampling [28] was used for the analysis of the mixed uncertainty quantification in the P-box. Each sample taken from epistemic distributions in the outer loop is only a possible value of the input with epistemic uncertainty, and it is supposed to be a fixed value in the inner loop. The simulation in the inner loop about a CDF of the system response quantities is considered as a conditional probability. When all of the epistemic samples are used, a series of CDFs have been completed.
In the past, the uncertainty quantification with the nested sampling were analyzed with the MCS/MCS in the P-box [25], [29]. However, in the presence of multiple uncertain parameters, especially when both epistemic and aleatory uncertainties need to be considered simultaneously, the computational cost to produce an accurate result is usually high. Hence, one of the motivations of this paper is to develop an alternative method, which should be effective and computationally efficient in quantifying both aleatory and epistemic uncertainties. Therefore, the NIPC method is employed for the inner loop and a two-phase MCS/NIPC method is introduced in the P-box for the flutter analysis with multiple uncertain variables in this paper.
Based on the characterization of uncertainty in system performance boundaries, the quantification of margins and uncertainties (QMU) have been used as one of the tools to facilitate analysis and communication of confidence for certification of complex systems [30]. An accurate QMU is crucial for the analysis of aeroelastic system which is highdemanding in performance requirement, and the uncertainty propagation process plays an important role in the procedures of the QMU. The concept of the QMU methodology was first proposed in 2003 for the safety assessment and decision-making process of nuclear weapons stockpile under the limited number of test data [31]. The main components such as performance gates, margins, and uncertainties [32], conceptual and computational basis [33] of the QMU methodology had been illustrated. Recently, there has been many discussions of the QMU application in many engineering problems such as radiative shock system [34], radioactive waste disposal [35], satellite power system evaluation [36], and missile reliability [37]. Thomas et al. [38] presented a practical QMU for complex spacecraft systems in which uncertainty exists in both the operating region and the performance requirement. Shah et al. [30] implemented Dempster-Shafer Theory of Evidence in the presence of mixed uncertainties in the reliability and performance assessment of complex engineering systems through the use of the QMU methodology. Xie et al. [39] proposed to use the QMU for simulation-based structural analysis of a pressure vessel with corrosion damage. But there is still little available work about using the QMU methodology in the stability analysis of aeroelastic system considering both the aleatory and epistemic uncertainties. Hence, another motivation for this paper is to construct a novel implementation framework of quantification of margins and mixed uncertainties (QMMU) for flutter speed analysis based on the results obtained by the MCS/NIPC method. After the introduction of the MCS/NIPC method and the QMMU, two simulation models, a twodimensional airfoil and a three-dimensional AGARD 445.6 wing model, are employed to demonstrate the effectiveness and efficiency of the proposed frameworks.
The remainder of the paper is structured as follows. Section II presents the concepts and theories including basics of flutter speed analysis and the two-phase MCS/NIPC method utilized for mixed uncertainty quantification. Section III introduces the QMMU framework. Detailed applications of the proposed approach to the 2D airfoil model and 3D AGARD 445.6 wing model are provided in Section IV. Finally, Section V concludes the paper by summarizing the current work.

A. BASICS OF FLUTTER ANALYSIS
Typical flutter equation of the undamped system is given as: where [M ] and [K ] represent the matrices of generalized mass and stiffness, respectively.
[A] is the aerodynamic influence coefficient matrix. V , q, and ρ are the flight speed, the generalized coordinates, and the atmospheric density, respectively. Generally, the unsteady aerodynamic matrix is a complex function in terms of Mach number and the reduced frequency k. The reduced frequency, k, is defined as: where ω is the cycle frequency (radian/s) and c is the reference chord length.
There are several methods used to solve the flutter equation, such as the V-g method, the p-k method, the k method, and the g method [40]- [42]. The first two methods mentioned above are adopted in this paper. In addition to analytical methods, the finite element software also widely used in the aeroelastic field. ZAERO [43], and MSC.NASTRAN [44] are recognized as the two most commonly used software for the flutter analysis, and the latter is chosen to conduct flutter analysis in this paper.

B. PROBABILITY BOX UNCERTAINTY QUANTIFICATION
P-box is an effective method for the mixed uncertainties propagation. A typical P-box is shown in Fig. 1, F(x) is a CDF from the simulation. Supposing F andF are non-decreasing functions in [0, 1], with F(x) ≤F(x) for all X ∈ R. Let [F,F] denotes the set of all non-decreasing functions F(x) in the interval [0, 1] written as F(x) ≤ F(x) ≤F(x). As can be seen from Fig.1, for a certain probability F(x j )(horizontal line), the right boundary F X (x) corresponds the upper valuex j of system response, and the left boundF X (x) corresponds the lower value x j . For a certain value x i (vertical line), F(x i ) is the lower probability whileF(x i ) is the upper probability. In the expression of a P-box, the likelihood of occurrence for system response quantity is not displayed as a precise probability, but an interval-valued probability. In other words, for a given value of the system response quantity, the corresponding probability can be expressed by an interval-valued probability. Similarly, for a given probability value, the corresponding predicted value of system response quantity is an interval value. For the case of mixed uncertainty quantification, the 95% confidence interval takes the maximum value of the upper 97.5% probability level and the minimum value of lower 2.5% probability level as the interval. The measurement is also illustrated in Fig. 1.
A two-phase MCS/MCS method is adopted in this study to get the P-box result [29]. The method generates exact values of uncertain variables randomly in the corresponding simulation to reach a result with proper accuracy by numerous iterations. The relationship between the number of samples (N s ) and the error of MCS (MCS e ) was given as [45]: The accuracy of the MCS increases as the number of samples increases. The number of samples is determined by the complexity of the relationship between input variables and output variables. Before using the MCS, different numbers of samples should be tested to determine the appropriate number of samples.

C. TWO-PHASE MCS/NIPC METHOD FOR MIXED UNCERTAINTY QUANTIFICATION 1) POINT COLLOCATION NON-INTRUSIVE POLYNOMIAL CHAOS
Polynomial chaos expansion is a powerful technique to provide a functional approximation of a computational model through its spectral representation on a suitably built basis of polynomial function [46]- [49]. Consider a random vector X ∈ R n with independent components, the polynomial chaos expansion of M (X) is defined as: where ψ i (X) is multivariate orthonormal polynomial and α i is the corresponding coefficient. For the i th -input random variable x i obeying Gaussian distribution: where ξ i is a function of the variable x i . E (x i ) and var (x i ) are the mean and variance of the variable x i , respectively. ξ = {ξ 1 , ξ 2 , . . . , ξ i , . . . , ξ n } is the N-dimensional standard random variable vector. Then the computational model is Theoretically, the number of terms of polynomial chaos expansion should be infinite [30], [50]. It is straightforward to define a ''standard truncation scheme'' given as: where p represents the order of the polynomial chaos and n represents the number of random dimensions. The proper order can be determined by comparison between the result of the NIPC method and the MCS method. Generally, the proper order is usually 2 or 3 for simple computational models and 4 to 6 for complicated models. Polynomial chaos expansions of Y is expressed as: Traditionally, the polynomial basis ψ i (ξ ) in (4) is built by using a set of univariate orthonormal polynomials φ (9) where i identifies the input variable with respect to which they are orthogonal, j and k represent the corresponding polynomial degree, f ξ i (ξ i ) is the i th input marginal distribution and δ jk is the Kronecker symbol. Note that the definition of the inner product can be interpreted as the expectation value of the product of the multiplicands.
Then the multivariate polynomial ψ i (ξ ) is assembled as the tensor product of its univariate counterparts: Due to the orthogonality relations in (9), it follows that the multivariate polynomials constructed are orthonormal: where δ αβ is an extension of Kronecker symbol. The classical families of univariate orthonormal polynomials and their corresponding distributions are given in Table 1. Detailed description of each classical families of polynomials are presented in [51].
There are several ways to calculate the coefficients α i of the polynomial chaos expansion for a given basis in (4). In general, an intrusive approach is required to modify the deterministic model and it is time-consuming and difficult when handling complex system [52]. To avoid the inconveniences in using the intrusive approach, non-intrusive approaches have been developed for the uncertainty propagation. There are two principal strategies that are used to calculate the polynomial chaos coefficients non-intrusively: the projection method [53] and the point-collocation NIPC [54]. The point collocation NIPC is briefly explained and implemented in this paper.
The point-collocation method first replaces the uncertain variables with their polynomial expansions given by (4). Then, N t vectors are selected randomly for a given polynomial chaos expansion with N t modes. The deterministic value is then calculated at the selected points as the left hand side of (4). Following this, a linear system of N t equations is formulated and then solved for the spectral modes of the random variables. This system is shown as: When the number of linear independent available samples is more than N t , the system is overdetermined and can be solved using the least squares approach. In the paper, the number of collocation points are determined twice as N t for point-collocation NIPC [55].  Fig. 2.
The following gives a brief outline of the procedure: 1. Set the number of samples for epistemic uncertainties as M ; 2. Set the number of samples for aleatory uncertainties as N ; 3. Determine M different values from epistemic uncertainties; 4. Choose the corresponding orthonormal polynomials, determine the proper order, select collocation points randomly, and construct the NIPC model; 5. Choose N random samples from aleatory uncertainties; 6. Evaluate the NIPC model to obtain the system response quantity by chosen samples; 7. Construct a CDF for the system response quantity using N samples; 8. Collect the M CDFs, and the result is an ensemble of CDFs, plotted as a P-box.
After the whole process, a series of cumulative distribution functions are plotted to give a P-box representation of the output variable. For the observed values, the largest value of the upper boundary and the smallest value of the lower boundary of the confidence interval are stored. The interval formed by the upper and lower boundaries is the confidence interval for the system response quantity involving aleatory and epistemic uncertainties. The confidence intervals with respect to different probability levels and the other important characteristics of the output variable can be obtained.

III. QUANTIFICATION OF MARGINS AND MIXED UNCERTAINTIES FRAMEWORK
Margins on system performance are important issue when uncertain variables exist. In a basic QMU framework, the QMU confidence ratio (CR) is defined as [56]:   where M and U represent the measure of the margin and uncertainty, respectively. If the confidence ratio CR is larger than 1, it means that the QMU metrics will be within the performance gate and indicates the performance boundary is safe when uncertainties are considered. Detailed applications of the QMU to various engineering fields differs from the definitions of the M and U , as showed in Fig. 3 [57].
In the representation of the QMU shown in (a) and (c), the operating region is allowed to have uncertainty (in the form of an interval). Margin is defined as the distance from the most conservative possible value in the operating region to ''the performance requirement''. In the representation of the QMU shown in (b), uncertainties exist on the operating region and the performance region. These are represented as two probabilistic distributions and margin is then defined as the distance between the means. In this paper, a practical QMMU framework is formulated and illustrated for the flutter speed boundary analysis with the aleatory and mixed uncertainties respectively, as shown in Fig. 4 and Fig. 5. Generally, the actual performance boundary is calculated by using a safety margin for the deterministic flutter speed. Within the definition of the QMMU in this paper, the uncertainty of the theoretical performance is considered, described by a Gaussian distribution.
For the aleatory uncertainty case, the results can be described by two curves as shown in Fig. 4, where P represents probability, VL and V f denote the actual and theoretical boundary of performance, respectively. β represents confidence interval, and β is 0.95 in this paper. U includes two parts, U 1 and U 2 , as shown in (14)- (16), while M is defined in (17).
For the mixed uncertainty case, this paper constructs a novel framework called QMMU to analyze the system stability, as shown in Fig. 5. V left and V right represent the left boundary and right boundary of the P-box result of the mixed uncertainty quantification. U includes two parts, U 1 and U 2 , as written in (18)- (20), and M is defined as (21). P represents probability, β represents the confidence interval, and β is also set to 0.95 in this paper.

IV. CASE STUDIES A. QMMU FOR THE TWO-DIMENSIONAL AIRFOIL 1) MODEL DESCRIPTION
A typical two-dimensional (2D) airfoil [10], as shown in Fig. 6, which has two degrees of freedom (pitching (α) and plunging (h)), is selected to demonstrate the superiority of the present methods for the flutter analysis with mixed uncertainties. G and E are the center of gravity and stiffness center,  b denotes the airfoil semi-chord, a is non-dimensional distance from mid-chord to elastic axis, x α is the non-dimensional distance from elastic axis to G, V represents the flight speed, k h and k α are heave spring-stiffness and torsion spring-stiffness, respectively. Experimental results show that the 2D airfoil has a low flutter speed [58]. Table 2 gives the parameter of the structure of the 2D airfoil and the flight condition for the flutter analysis, where ρ is the density of air, m denotes the mass of 2D airfoil, I α is the rotational inertia [59]. Flutter analysis is conducted by using the p-k method and the flutter speed V f calculated from the deterministic system is 39.48 m/s.

2) QUANTIFICATION FOR MULTIPLE ALEATORY UNCERTAIN VARIABLES
In this section, the aleatory uncertain variables the stiffness parameters k α and k h are considered while the epistemic uncertainty is neglected. Here, the two parameters are modeled as independent Gaussian random variables, and the variations of the random variables are determined by a coefficient of variation (COV) of 5%. They are shown in Table 3.   The multidimensional Hermite polynomials are used in the NIPC calculations.
The effects of the aleatory uncertainty on the flutter speed are quantified by using the MCS (10000 simulations) and second-order Point-Collocation NIPC (only 12 points needed according to (7)). Obviously, the CDF and PDF of the NIPC model with the polynomial degree of 2 agree well with the ones obtained by MCS as shown in Fig. 7 and Fig. 8. Table 4 shows the mean, standard deviation, 95%confidence interval and the computational time associated with the two methods. The first three results can also be observed visually from the CDF plots and PDF plots. Moreover, as can be seen from the table, the computational time of the MCS is approximately 10 times that of the Point Collocation NIPC.

3) QUANTIFICATION FOR MIXED UNCERTAIN VARIABLES
In engineering applications, aleatory uncertainty and epistemic uncertainty exist in one system. Consequently, a more comprehensive flutter speed boundary prediction is desired when taking both uncertainties into account simultaneously.
Both aleatory and epistemic uncertainties are considered in this section. The parameters k α and k h are still treated VOLUME 8, 2020   as aleatory uncertain variables while the mass m and the rotary inertia I α are treated as epistemic uncertain variables, because the accurate PDFs of these two parameters cannot be obtained, nor the accurate mean values and variances. The parameter m is a value between [11.15, 13.63], and the parameter I α is a value between [0.059, 0.072]. All of the uncertain variables are summarized in Table 5.
There are 100 × 10000 iterations that have been carried out in this section. The final results are shown in Fig. 9, and the flutter speed is in the form of a P-box, which could represent  the output caused by both the epistemic and aleatory uncertainties. It is notable that only 57.58 seconds are taken for the two-phase MCS/NIPC (second-order) method to generate the final results. The computational efficiency is significantly higher than that of using the two-phase MCS/MCS method (5078 seconds). When the computational model is more complicated, the MCS/MCS method needs numerous simulations which is too time-consuming to use, but the MCS/NIPC is an alternative method, which is effective and computationally efficient.
Each iteration in the outer loop generates a CDF based on the aleatory uncertainty analysis in the inner loop. Consequently, the width between the leftmost CDF and rightmost CDF reflects the contribution of epistemic uncertainty. The 95% confidence interval, which is more concerned in engineering process, is obtained by taking the maximum value of the 97.5% probability level and the minimum value of the 2.5% probability level. For this problem, the 95% confidence interval is given as [36.94, 42.06]. The probability levels and the corresponding probability intervals for the P-box are illustrated in Table 6. Moreover, for the exact flutter speed 39.48 m/s which is generated by deterministic computation, the corresponding probability level is [0.37, 0.58], obtained by taking the maximum value of the far-left CDF and the minimum value of the far-right CDF.

4) QMMU FOR FLUTTER EVALUATION OF 2D AIRFOIL
First, the QMU application for aleatory uncertainty quantification of 2D airfoil is performed. From previous sections, the distribution of flutter speed under uncertain conditions,  obtained by means of modeling, simulation and uncertainty quantification, is the theoretical flutter speed boundary, as shown in Fig. 7. Considering there are uncertainties such as calculation error, the safety margin is 15% and the variance is determined by a coefficient of variation of 2.5% for the deterministic (or actual) flutter speed. The actual flutter speed is assumed to satisfy the normal distribution as shown in Fig. 10. The QMU metrics can be obtained by applying the framework and equations in Section III. From the results, as shown in the Table 7, the confidence ratio for aleatory uncertainty case CR 1 is smaller than 1 in terms of these uncertainties, which means that the ''Margin'' is not sufficient to cover the ''Uncertainties''. As a result, the flutter speed boundary with 15% safety margin is not safe enough.
Then the QMMU analysis for mixed uncertainty quantification of 2D airfoil is conducted. Similarly, the theoretical flutter speed boundaries are the same as the results in Fig. 9.The actual flutter speed, which is assumed to satisfy normal distribution with the safety margin of 15% and the coefficient of variation of 2.5%, is shown in Fig. 11. By applying the framework and equations in Section III, the confidence ratio for mixed uncertainty case CR 2 is 0.43, which means ''Margin'' is smaller than the ''Uncertainties''. All the results of the QMMU metrics are summarized in Table 8. Obviously, CR 2 is smaller than CR 1 due to the consideration of epistemic uncertainty. Both for the aleatory uncertainty case and mixed uncertainty case, the values of CR are smaller than 1, indicating that a redesign of flutter speed limits maybe required to make the system more reliable.  To further illustrate the application of the MCS/NIPC method and the QMMU framework to the flutter analysis with mixed uncertainties, the AGARD 445.6 wing model, a swept-back wing at transonic speed, is selected in this section. The quarter-chord sweep angle is 45 degrees, the aspect ratio is 1.65, the taper ratio is 0.66, the wing root chord is 558.8 millimeters, and the wing semi-span is 762 millimeters with a symmetric NACA 65A004 airfoil section. The geometry of the wing is shown in Fig. 12 [60].
The deterministic flutter analysis of the AGARD 445.6 wing is computed by using the flutter solution sequence, SOL 145 in MSC Nastran at Mach number 0.5. Table 9 gives the material parameters of the model. E 11 , E 22 represent elasticity modulus in material principal direction, x and y-axial direction, G 12 denotes the shear modulus, υ is the Poisson's ratio, and ρ stands for the density of the model, respectively. r, which represents the ratio of the air density in the current flight altitude to the one in a certain altitude, is 0.3486. The deterministic flutter speed V f is 178.65 m/s solved by p-k method.

2) QUANTIFICATION FOR MULTIPLE ALEATORY UNCERTAIN VARIABLES
The above material parameters could not be accurate in the complex flight environment due to the variations of temperature, pressure, flight altitude and so on. The point collocation NIPC method is adopted to address the aleatory uncertainty propagation in the AGARD 445.6 wing model. Here, conventionally the material properties E 11 , E 22 , G 12 are treated VOLUME 8, 2020    as aleatory uncertain variables and Gaussian distribution is chosen to model the uncertain parameters. Therefore, multidimensional Hermite polynomials are used according to Table 1 in the NIPC calculation. The three parameters are modeled as independent Gaussian random variables, and the variations of those variables are determined by a coefficient of variation (COV) of 5%, as shown in Table 10.
To develop a functional approximation of the computational model, the order of the polynomial chaos expansions is decided with respect to the number of collocation points as given in Table 11. The sampling method used here is Latin Hypercube sampling (LHS).
The results obtained by using the NIPC method with different orders in the polynomials are compared with the result obtained by the MCS method (1000 simulations). Obviously, as is shown Fig. 13 in the curves of the first-order NIPC  and the second-order NIPC do not agree well with the curve of the MCS result. However, there is not much difference between the result of the third-order NIPC and the MCS, which leads to a conclusion that the NIPC with the third-order polynomials can describe the propagation of the aleatory uncertainty for this model as an alternative approach. As higher-order polynomials increase the computational cost, third-order NIPC is chosen in this paper for model simulations. Table 12 summarizes the specific characteristics of flutter speed distributions of four curves in Fig. 13, revealing how the aleatory uncertainty in the material properties influences the flutter speed boundary. As seen from the result obtained using the third-order NIPC method, 95% confidence interval of the flutter speed is [172.93 185.60] when taking the three uncertain material parameters into account. It is noted that the coefficient of variation of the outputs generated by any of four methods is less than that of the inputs.

3) QUANTIFICATION FOR MIXED UNCERTAIN VARIABLES
Both the material parameters and flight condition parameters would change in the practical applications. Here, the parameter r, which is used to represent the ratio of the air density in the current flight altitude to the one in a certain altitude,  can be considered as an epistemic uncertain variable due to the variation in the flight altitude. Generally, the uncertain variable that cannot be described by an exact distribution due to the lack of knowledge can be regarded as interval value as well. Here, the material density ρ is treated as epistemic uncertain input variable in a form of interval. Therefore, there are 5 uncertain input variables and they are summarized in the Table 13.
In this part, 100 × 1000 iterations of numerical simulation have been conducted by using the two-phase MCS/NIPC (third-order) method. The P-box presentation of flutter speed is given in Fig. 14, in which the impact of the epistemic uncertain variables is illustrated by the CDFs, and the impact of the aleatory uncertain variables is illustrated by the width of each curve between the far-left side and the far-right side. It is clearly seen that the flutter speed has a nonignorable fluctuation in the case of the mixed uncertainty quantification which contains 3 aleatory uncertain variables  and 2 epistemic uncertain variables. For the whole P-box presentation, the 95% confidence interval is [165. 51,198.09] by taking the maximum value of the 97.5% probability level and the minimum value of the 2.5% probability level. The interval reveals the joint contribution shared by both aleatory and epistemic uncertainty. Then, the probability levels and corresponding probability intervals for the P-box are illustrated in Table 14, illustrating how the epistemic uncertainty influences the flutter speed. For the left-most CDF of the P-box, the 95% confidence interval is [165. 51,176.84] and for the right-most CDF is [181.91, 198.09]. The two intervals depict how aleatory uncertainty in the material parameters influences the flutter speed when the material density and flight conditions are given. All of the results demonstrate that both the aleatory uncertainty and epistemic uncertainty have a significant effect on the flutter speed boundary.

4) QMMU FOR FLUTTER EVALUATION OF AGARD 445.6 WING MODEL
Similar to the previous two-dimensional airfoil problem, the next step is to perform the QMU framework for aleatory uncertainty quantification of AGARD 445.6 wing. As well as the safety margin of 15% and the variance determined by a coefficient of variation of 2.5%, the actual flutter speed is assumed to satisfy the normal distribution with the CDF of actual boundary shown as dashed line in Fig. 15. The theoretical boundary is the same as the third-order NIPC one in Fig. 13. Using the framework and equations in Section III, the confidence ratio for aleatory uncertainty case CR 1 is 1. 35. It means that the ''Margin'' is greater than the VOLUME 8, 2020   ''Uncertainties'' for the case where only aleatory uncertainty exists and the flutter speed boundary is safe. All the results of the QMU metrics are summarized in Table 15.
For the QMMU application for the mixed uncertainty quantification of AGARD 445.6 wing as shown in Fig. 16, the theoretical flutter speed boundaries are consistent with the results in Fig. 13. By applying the framework in Section III, the confidence ratio for mixed uncertainty case CR 2 is 0.22. This means the ''Margin'' is smaller than the ''Uncertainties'', indicating that a redesign of flutter speed limits maybe required. All results of the QMMU metrics are summarized in Table 16. Obviously, CR 2 is much smaller than CR 1 because of the additional consideration of the epistemic uncertainty. As CR 1 is larger than 1 and CR 2 is smaller than 1, it implies that considering the epistemic uncertainty of engineering systems is necessary.

V. CONCLUSION
The present paper proposes an efficient mixed uncertainty quantification method named the two-phase MCS/NIPC method and introduces a practical QMMU implementation framework in the presence of mixed uncertainties for flutter speed boundary computation. To reduce the computational cost without a loss of accuracy in the predictive results, the point collocation NIPC is employed for the uncertainty quantification. The methods are demonstrated by computing the flutter speed boundary for a typical two-dimensional airfoil and a three-dimensional AGARD 445.6 transonic wing model.
The proposed two-phase MCS/NIPC method is practical and efficient for the flutter analysis when the systems suffer from the multiple aleatory and mixed uncertainties due to the variation of material parameters and flight conditions. Aleatory uncertainty quantification and mixed uncertainty quantification are considered, respectively. From the results of aleatory uncertainty quantification, the computational efficiency of the NIPC method could be demonstrated. The results of mixed uncertainty quantifications described by a P-box demonstrate the applicability of two-phase MCS/NIPC method in the flutter speed boundary prediction for complex models.
The comparisons of the QMU and QMMU applications to the two models are presented to illustrate the significant influence of epistemic uncertainty, indicating that mixed uncertainty quantification for aeroelastic systems are necessary. The QMMU procedure could give some help to determine the safety margin as it can give an analytical conclusion of system performance based on QMMU metrics, showing the potential for the evaluation of stability of complex engineering systems. There are many other methods for P-box estimation which could be used in aeroelastic systems, and they will be considered in follow-up research.