Bifurcation Analysis with Aerodynamic-Structure Uncertainties by the Nonintrusive PCE Method

An aeroelastic model for airfoil with a third-order stiffness in both pitch and plunge degree of freedom (DOF) and the modified Leishman–Beddoes (LB) model were built and validated. The nonintrusive polynomial chaos expansion (PCE) based on tensor product is applied to quantify the uncertainty of aerodynamic and structure parameters on the aerodynamic force and aeroelastic behavior. The uncertain limit cycle oscillation (LCO) and bifurcation are simulated in the time domain with the stochastic PCE method. Bifurcation diagrams with uncertainties were quantified. The Monte Carlo simulation (MCS) is also applied for comparison. From the current work, it can be concluded that the nonintrusive polynomial chaos expansion can give an acceptable accuracy and have a much higher calculation efficiency than MCS. For aerodynamic model, uncertainties of aerodynamic parameters affect the aerodynamic force significantly at the stage from separation to stall at upstroke and at the stage from stall to reattach at return. For aeroelastic model, both uncertainties of aerodynamic parameters and structure parameters impact bifurcation position. Structure uncertainty of parameters is more sensitive for bifurcation. When the nonlinear stall flutter and bifurcation are concerned, more attention should be paid to the separation process of aerodynamics and parameters about pitch DOF in structure.


Introduction
Hopf bifurcations can occur in aeroelastic models that are nonlinear just in the structural stiffness operator (e.g., panel LCO), but, nonlinear aerodynamics is, in most cases, the critical component of computational aeroelasticity that must be enhanced to enable dependable predictions of aeroelastic response and variability. Dynamic stall which is due to nonlinear phenomenon for aerodynamics is one of the most important reasons that leads the change in bifurcation. Calculating nonlinear dynamics correctly is the key to the investigation of LCO and bifurcation. For nonlinear dynamics, CFD method is a direct way to calculate which have a good accuracy. However, when the problem is complex and many conditions need to be calculated, using CFD method may cost much time [1]. There are other methods to calculate nonlinear aerodynamics such as Leishman-Beddoes (LB) [2]. This method can improve the calculation efficiency greatly. It is a semiempirical model and many of the model parameters are selected subjectively; there may be some variations of these parameters. These parameter uncertainties in the aerodynamic model may lead to conservative or optimistic prediction of stall flutter [3,4].
Since uncertainty always exists in both aerodynamics part and structure part when a theoretical method is applied to the aeroelastic stability analysis, aeroelastic uncertainty quantification (AUQ) is somewhat unavoidable in the theoretical analysis [5]. Dai and Yang Reviewed the methods applied to aeroelastic uncertainty quantification [6], which is focused on investigating the effect of parameter uncertainty on the aeroelastic stability property [7,8]. The methods can be divided into robust one [9] and probabilistic one [10]. When the uncertainty is a probabilistic variable, the stochastic aeroelastic analysis should be conducted to obtain both the stability boundary and its distribution [8,11,12]. The nonintrusive polynomial chaos expansions are to determine the evolution of uncertainty in a dynamic system, when a probabilistic uncertainty is embedded in the aeroelastic system [13,14]. Beran et al. introduced this method to the nonlinear aeroelastic analysis, which has been proven to be a success [15]. 2 International Journal of Aerospace Engineering Some related work has been carried out, including uncertainty quantification in aerodynamics and stochastic basis construction, Badcock et al. [16][17][18]. In AUQ, the parameter uncertainty exists both in the structural dynamics model and the aerodynamics [19,20]. Most of these works are focused on the structural uncertainty in structure [21,22]. Little work is concerned with aerodynamic uncertainty and structure uncertainty interactions in aeroelasticity [23,24].
This current work is focused on the aerodynamic uncertainty quantification of aeroelastic stability containing limit cycle oscillation analysis and flutter analysis, subject to stochastic LB aerodynamic parameters at low Mach number. The modified unsteady aerodynamic model at low Mach number will be established, and uncertainty of parameters for aerodynamic model will be taken into consideration. Then the stochastic polynomial chaos expansion method for stall flutter will be investigated.

Deterministic Aerodynamic
Model. The Leishman-Beddoes (LB) dynamic stall model is a popular model that has been used to investigate dynamic stall of aerodynamics. The standard LB model is applicable above the Mach number of 0.3. But, it is not suitable in the incompressible flow. However, a huge effort has been undertaken for the aeroelastic characterization of next-generation aircraft; particularly, the low Mach number is important for high altitude very-long endurance solar power HALE UAVs [25,26]. What is more, in most of stall flutter application cases, such as wind tunnel state, the Mach number is lower than 0.3. Hence, in this paper, a modified dynamic stall model for low Mach numbers model was built based on the work of Leishman and Beddoes and Sheng et al. [2,27]. Under the conditions of low Mach number, an additional time lag is required for the disturbed flow, during which the disturbed flow can develop into vortex that is strong enough to cause the dynamic stall. Taking this effect into consideration, the additional lagged value , which is the substitute value of , is introduced.
where is an intermediate variable, is the time lag constant, and is the distance traveled by airfoil in semichords.

Modification of Normal Force and Pitching Moment
Coefficients at Upstroke. Vortex forms and detaches near the airfoil leading edge, and the flow of the separated shear layer still attaches to the upper surface at low Mach number. Hence, it results in additional overshoot in normal force at low Mach numbers. The overshoot for normal force coefficient Δ V is given as follows: where 1 is a coefficient related to airfoils; is separation location in terms of chord; is delayed separation location of in the original LB model; represents the shape function of normal force due to vortex, which is given by where is nondimensional time; V is time constant of vortex traveling over chord; V is vortex passage time constant. The additional pitching moment coefficient Δ V is also the effect of vortex convection over the airfoil. It is as follows: where 2 is the coefficient dependent on airfoil. V is nondimensional time during vortex passage.

Modification of Normal Force and Pitching Moment
Coefficients during Return. During return, there is also an overshoot for normal force coefficient. A value of min 0 was given. It is assumed that it is the start of the convection process of overshoot. A value of min was given to define the end of the convection process of overshoot. Their relationship is as follows: where is delay constant for reattachment process and is reduced pitch rate. The overshoot for normal force coefficient at return process is given as follows: where represents the shape function of normal force due to vortex at return process, which is given by where is a nondimensional time variable at return process.

Deterministic LCO Analysis in the Time Domain.
The aeroelastic system considered here is a wing section with pitch and plunge motions, which is shown in Figure 1.
Considering the unsteady aerodynamic force, the motion equation of the wing section is written as where is the airfoil pitch angle. ℎ is the vertical displacement. is the mass per unit length. ℎ and are spring stiffness in bending and torsion. ℎ and are damping in bending and torsion. is static mass moment. is the polar moment of inertia about 1/4 chord. ℎ and are external aerodynamic force and moment, respectively, which will be calculated by the above modified LB model at low Mach number. There expressions can be given as follows: where and are lift coefficient and moment coefficient, respectively, is the chord length, is the span length, ℎ is the distance between 1/4 chord and elastic axis, is air density, and is the free stream velocity.
The angle of attack and pitch ratio which are used in LB model can be given as follows: ] .
For a nonlinear aeroelastic system, simulation in the time domain is a common and good approach to investigate its time response, such as the limit cycle oscillation phenomenon. In the current study, a fixed step 4th-order Runge-Kutta algorithm is applied, both for the deterministic solution and in the uncertain cases.

Nonintrusive Polynomial Chaos Method for Nonlinear
Aeroelastic System. The LB model has many aerodynamic parameters; these aerodynamic parameters are different at different Mach numbers and airfoils. All the parameters are gotten from experiments, which are the ensemble average of some 50 pitch cycles [28]. Hence, stochastic uncertainties widely exist in parameters, which are gotten from experiments. In modified LB model, the number of parameters is more than original model. The uncertainty of these parameters brings uncertainty to the unsteady aerodynamics which will add uncertainty to airfoil aeroelastic system. Specifying uncertainty of aerodynamics in airfoil aeroelastic system, 16 parameters in modified low Mach numbers LB model were investigated. They were the 14 parameters in original LB model: , 1 , 1 , 2 , 0 , 1 , 2 , 0 , , 1 , , , V , and V and the two parameters in modified low Mach number LB model: and . All the uncertainties are noted as = [ 1 , 2 , . . . , ]; here = 16. Then responses in the time history are written as x( , ), which will be also a stochastic vector.
PCE estimates its coefficients using the approaches random sampling or tensor-product quadrature [29]. To build the uncertainty model, a tensor product of one-dimensional rules was used. Since different distribution of parameter error may affect the result of calculation, here it is assumed that the distribution of parameter error was normal or uniform which was calculated separately.
To construct the stochastic orthogonal basis of outputs, the expansion of the aeroelastic response can be expressed: where ( ) are multidimensional orthogonal polynomial basis which are stochastic part. The expansion is truncated to term which is determined by the number of variables and the order of polynomial ( ).
This is referred to as total-order expansion. In order to compute the stochastic output equation, it is required to solve a problem with unknowns and equations. Therefore, it is required to run the model for times in order to evaluate the deterministic functions which are used to evaluate the mean and variance 2 . The statistics of the random solution are given by With ⟨ ⟩ denoting the inner product It is possible to evaluate the deterministic functionsx ( ) ( = 0, 1, . . . , − 1) by using the orthogonality of ( ).
where each inner product involves a multidimensional integral over the support range of the weighting function. Full tensor product is to employ a tensor product of one-dimensional quadrature rules and was used to estimate coefficient. In the tensor-product case, Gaussian abscissas were chosen, the zeros of polynomials that are orthogonal with respect to a density function weighting. For uniform and normal distribution, Gauss-Legendre and Gauss-Hermite were used separately.
One-dimensional quadrature rules are as follows: where 1 is the number of quadrature point , for onedimensional quadrature rules. , is the weight of quadrature point. Then, for multidimensional quadrature rules, the equation changes as follows: For calculating uniform distribution, the Legendre polynomials are applied. In this polynomial series, the range of random variables is located in [−1, 1]. The probability density function of each variable is 1/2 in the range of [−1, 1].
The Legendre polynomial bases are written as Le 0 = 1, x ( ) are deterministic functions to be evaluated and Le ( ) are stochastic part tabulated for uniform distribution. Then, (15) can be written aŝ Then, (17) can be used to calculate the coefficient. The information of quadrature point and its weight for Gauss-Legendre used in this paper is shown in Table 1.
is the number of quadrature points, which can be gotten by using = + 1.
Then, the aeroelastic response of the stochastic process can be calculated directly according to (11). The numerical sampling of random variables is still needed in the above nonintrusive polynomial chaos method. It is advantageous compared with the standard Monte Carlo simulation. Considering the stochastic variation for the aerodynamic parameters, they are represented by the Wiener expansion where 1 , 2 , . . . , are the random variables satisfying a uniform distribution and 1 , 2 , . . . , are their corresponding uncertainty bounds, respectively. Consequently, the response x in (11) can be rewritten as In (21), 0 represents the mean value of the response and is the corresponding fluctuation of the random components. Using the first-order representation of Legendre polynomials, the response can be written as 16 16 .
For calculating normal distribution, the Hermit polynomials are applied. In this polynomial series, the range of random variables is located in [−∞, +∞]. The probability density function of each variable is The Legendre polynomial basis is written as International Journal of Aerospace Engineering 5 Then (15) becomes as follows: The information of quadrature point and its weight for Gauss-Hermite which is used for normal distribution is shown in Table 2.

Deterministic Model Validation
To make sure the model has a good accuracy, which can be used to investigate the uncertainty of airfoil aeroelastic system, validation is conducted for modified low Mach number model and airfoil aeroelastic system, respectively, in this section.

LB Model Validation at Low Mach Numbers.
To validate the modified LB model at low Mach number, the aerodynamic coefficients of NACA0012 were investigated. The results of this paper are compared with results of experiment, Sheng, and original model [27]. The freestream velocity is 0.1 Mach. The results at the condition of deep dynamic stall (angle of attack is = 15 ∘ + 10 ∘ sin ) are compared. The compared results are shown in Figures 2 and 3.
From Figures 2 and 3, it can be seen that the modified low Mach number LB model in the current work can match results of experiment mostly. The modified low Mach number LB model can calculate overshoot at upstroke and return correctly. As a deterministic aerodynamic calculation model, the modified low Mach number LB model has a good accuracy, but there still are some test data which are not covered by the deterministic aerodynamic model. To improve the accuracy of calculation, the uncertainty of aerodynamics needs to be taken into consideration.

Validation for LCO Analysis.
Reference [30] presents a set of experiments conducted on NACA0012 airfoil undergoing the stall oscillations at low Mach number, and they are chosen as the numerical example to verify and to investigate the deterministic airfoil aeroelastic system model mentioned previously. The important experiment parameters are shown in Table 3.
To validate the deterministic model of aeroelastic system, the pitch bifurcation diagram and the LCO trajectories were calculated, since the aeroelastic system studied here is   model is an expression of linear stiffness. A third-order stiffness was added in both pitch and plunge degree of freedom (DOF). For present modified model, the restoring moment associated with the torsional spring was expressed as where 3 is the dimensionless parameter governing the third-order term. In DOF of plunge, similar expression of force can be given as The result of present work was compared with the result of experiment and the result of Song [31], which are shown in Figures 4 and 5.
From Figures 4 and 5, it can be seen that the deterministic results for modified model in the present work are well coincided with the experiment results which fits better than original model in amplitude. The results of phase plane between anḋfor both modified model and original model fit much better to the experiment than the result of Song. However, there are still some test data which are not covered by the result data of deterministic model. Hence, for improving the accuracy, the uncertainty in airfoil aeroelastic system should be investigated. This work can be carried out by using the same deterministic LCO model in the current work. It is therefore called the nonintrusive method.

Numerical Examples for Uncertainty Quantification
In current work the uncertainty of unsteady aerodynamics was investigated first. The impact of uniform and normal distribution for parameters on the calculation result were investigated separately and compared together. Then both unsteady aerodynamic and structure uncertainty were taken into consideration to investigate the dynamic stall flutter at low Mach number.

Uncertainty Quantification of Unsteady Aerodynamic
Force. In this section, for uniform distribution of parameters, the first-order Legendre polynomial basis was used to build nonintrusive PCE model for unsteady aerodynamics. For normal distribution of parameters, the first-order Hermit polynomial basis was used. For normal distribution, it is assumed that the distribution is (0, (0.0333) 2 ) and (0, (0.0666) 2 ), the mean = 0, and standard deviation is = 0.0333 and = 0.0666. The corresponding uncertainty bounds of aerodynamic parameters for uniform distribution are set as 1 = 0.1 and 1 = 0.2, and they are calculated, respectively. The plots for two kinds of distribution are compared and shown in Figure 6.
The result for Monte Carlo method was used as a reference method to compare with the result of nonintrusive PCE method for two kinds of distribution in the investigation of uncertainty for unsteady aerodynamics. Figures 7 and 8 show the mean of normal force coefficient and moment coefficient for two kinds of distribution. For uniform distribution, corresponding uncertainty bounds for Monte Carlo and nonintrusive PCE method are 0.2. For normal distribution, the standard deviation is 0.0666 for Monte Carlo and nonintrusive PCE method. Figure 9 shows the variance of normal force coefficient ( ) at upstroke and at return for two kinds of distribution. From Figures 7 and 8, it can be seen that when the uncertainty was added in the model, the mean results of Monte Carlo method with 30000 times simulation and nonintrusive PCE simulation will be changed. The results under two kinds of distribution are different. Under two kinds of distribution, the results of nonintrusive PCE method match the results of Monte Carlo method well, which can illustrate that nonintrusive PCE method in current work is correct. From Figures 7 and 8, it also can be seen that two kinds of  distribution have almost unanimous uncertainty bond which can be seen from Figure 6(b); since normal distribution has a smaller deviation ( = 0.0666) than uniform distribution, the mean of normal distribution is closer to deterministic than uniform distribution. It can be seen that the uncertainty impacts the stage of after separation and the stage of reattach obviously. The mean of Monte Carlo method and nonintrusive PCE method fit each other well even at overshoot during upstroke and return for two kinds of distribution. However, it can be seen from Figure 9 that the variance of normal force coefficient for two kinds of distribution especially at overshoot during upstroke is different. The number of variables in the calculation of uncertainty for aerodynamics is 16 which leads the method of PCE to not have a higher calculating efficiency than method of MCS. Figure 10 shows the normal force coefficient bound by 30000 times of MCS and nonintrusive PCE times for uniform distribution with corresponding uncertainty bound 0.1 and normal distribution with standard deviation 0.0333. From the figure, it can be observed that the uncertainty for parameters can affect the part of stall mostly. The normal force coefficient bound near stall during upstroke and reattach during return will expand obviously. However, the bound at stage of attach changes slightly. This is mainly because most parameters with uncertainty are related to the stage of stall.
The bound of the normal force coefficient due to normal distribution is smaller than uniform distribution at the stage of stall and reattach. The test data can be covered mostly by the bound when the distribution is uniform distribution.
By using the nonintrusive PCE method, surrogate expression of normal force coefficient and moment coefficient for two kinds of distribution in (23) was gotten. Parameters were nondimensionalized by and separately. Figures 11 and 12 show uncertainty parameters in normal force coefficient for two kinds of distribution and Figure 13 shows uncertainty parameters in moment coefficient for uniform distribution. In the three figures mentioned above, parameters with obvious sensitivity were labeled with wide lines. It can be observed from Figures 11 and 12 that the sensitivity of parameters is different during different stage; for two kinds of distribution they have similar tendency for different stage but the amplitude is different which can be translated by    19) and (25). and 1 are sensitive for both upstroke and return, especially at the return stage. It also can be seen that parameters with larger sensitivity values will affect stage of stall at upstroke. , V , and V are more sensitive parameters for during upstroke than others, while, during return, and become more sensitive than other parameters. From Figure 13, it can be observed that parameters which are sensitive for are also sensitive for . However, the number of parameters which are sensitive for is more than those for , such as 0 , 1 , and 2 . This is consistent with LB model. Since stall flutter is largely affected by dynamic stall and these parameters are related to dynamic stall, it is very important to get the relationship between these parameters and the force and moment.
From the parameters shown above, it can be seen that the parameters which are sensitive to also can effect obviously. For decreasing the calculating load, four parameters which have obvious effect on were filtered to use for investigating the stochastic uncertainty of LCO. The four parameters are and 1 , , and which will contribute most to normal for coefficient and moment coefficient . For two kinds of distribution, the boundary result of normal distribution is more conservative than uniform distribution. The tendency for sensitive parameters in two kinds of distribution International Journal of Aerospace Engineering is similar. Based on the reasons above, when investigating uncertainty of LCO and bifurcation, uniform distribution of parameters was used in calculation.

Uncertainty Quantification of LCO and Bifurcation.
In order to investigate the effect of uncertainty for bifurcation due to aerodynamics, the LCO needs to be calculated first.
To have an exact judgment for LCO, enough calculation cycles are needed. In the current work, when calculating the response of the pitch angle, the time step is 0.001 and the step number is 60000. Since uncertainty variables in calculating LCO are just 4, the calculation efficiency of method of PCE is higher than MCS method. For PCE method in current work, it just uses 16 times of simulation to calculate bifurcation diagram which cost 7053 s; however, for MCS method, an acceptable result needs 5000 times of simulation which cost 2204062 s which is not acceptable. Since stall flutter is a typical nonlinear aeroelastic phenomenon, not only aerodynamic parameters but also structure parameters have effect on it. Only taking uncertainty of aerodynamic parameters into consideration is insufficient. So the model with four aerodynamic parameters mentioned above was marked as model "Aerodynamic." Another model with four structure uncertainty parameters was marked as model "Structure." In model "Structure," four uncertainty parameters are torsional spring stiffness and damping and spring stiffness ℎ and damping ℎ . Figures 14, 15, and 16 show the response for the amplitude of the pitch angle at three different velocities which were calculated by the deterministic aeroelastic model. In Figure 14(a), the response of pitch angle converges with time. Figure 15(a) shows a standard LCO which occurs at the speed of 12.5 m/s when time passes 35 s. When the flow velocity continues to increase, the amplitude of LCO will reach a new value. Here it can be seen from Figure 16(a) that the LCO amplitude at = 16m/s is larger than the one at = 12.5 m/s. Using the judgment rule defined above, bifurcation diagram with stochastic uncertainty was gotten. Figure 17(a) shows the mean value for stochastic uncertainty bifurcation diagram of pitch angle amplitude for 1 = 0.1 and 1 = 0.2 in model "Aerodynamic" which are compared with the deterministic value and Figure 17(b) shows sensitivity of parameters for 1 = 0.2 and probability density function (PDF) of bifurcation. In order to explain the problem clearly, the parameters were nondimensionalized by divided peak in the bifurcation diagram. It can be seen from the figure that stochastic uncertainty of aerodynamics will change the mean value of bifurcation velocity. The mean bifurcation velocity for 1 = 0.2 is = 10.8 m/s which is smaller than the mean value of 11.5 m/s with the corresponding uncertainty bound = 0.1. The stochastic uncertainty of aerodynamics mainly affects the bifurcation point. In the four parameters chosen form aerodynamic model, only three of them are sensitive parameters for bifurcation. For parameters, the most sensitive stage is PDF from 0 to 1. The three sensitive parameters for bifurcation in model "Aerodynamic" are chosen as parts of parameters to build final model which contains uncertainty of both aerodynamic and structure. Figure 18(a) shows the mean value for stochastic uncertainty bifurcation diagram of pitch angle amplitude for 1 = 0.1 and 1 = 0.2 in model "Structure" which are compared with the deterministic value and Figure 18(b) shows sensitivities of parameters for 1 = 0.2 and PDF of bifurcation. The parameters in Figure 18(b) were also nondimensionalized. From Figure 18(a), uncertainties of structure parameters also can change the bifurcation position. The mean bifurcation velocity for 1 = 0.2 is = 10.8 m/s which is smaller than the mean value of 11.3 m/s with the corresponding uncertainty bound = 0.1. In Figure 18(b), it can be seen that, in the four uncertainties of parameters, two parameters about pitch are sensitive for bifurcation and pitch amplitude. Then these two parameters and were chosen as uncertainties parameters of structure to take part in building model "Aerodynamic + Structure" which contains both aerodynamic uncertainty and structure uncertainty.
Using the three uncertainties parameters chosen from model "Aerodynamic" and two parameters chosen form model "Structure", a model contains both aerodynamic uncertainty and structure uncertainty was built which was marked as "Aerodynamic + Structure" to investigate comprehensive impact for both aerodynamic and structure. Figure 19 shows the bifurcation diagram of pitch angle amplitude in "Aerodynamic + Structure"; it can be seen that due to the collective effect of aerodynamic uncertainty and structure uncertainty, the bifurcation velocity is smaller than model with single aerodynamic uncertainty or structure uncertainty whether the uncertainty bound is 0.1 or 0.2.
Statistics information of bifurcation velocity for three models was shown in Table 4. From the table it can be seen that for bifurcation position both aerodynamic uncertainty and structure uncertainty are very important. In "Structure" model, PDF from 0 to 1, the velocity covers a lager range than "Aerodynamic" model. Because of the collective effect 10 International Journal of Aerospace Engineering     for aerodynamic uncertainty and structure uncertainty, bifurcation occurs in advance. For models "Aerodynamic" and "Structure," it can be seen that bifurcation bound variation has a smaller range than the uncertainty bound: when 1 = 0.1, the uncertainty bound of parameters is 0.2, and the bound variation is 0.05 in "Aerodynamic" model and 0.083 in "Structure" model; when 1 = 0.2, the uncertainty bound of parameters is 0.4, and the bound variation is 0.108 in "Aerodynamic" model and 0.158 in "Structure" model. However in "Aerodynamic + Structure" model, when 1 = 0.1, the bound variation is 0.108; when 1 = 0.2, the bound variation is 0.241. Figure 20 is the sensitivities of parameters in "Aerodynamic + Structure"; it can be seen that structure uncertainty parameter is the most sensitive parameter for bifurcation; then aerodynamic uncertainty parameter which is delay constant for the flow separation point due to unsteady effect is also a big sensitive parameter second only to . The delay constant of the separation point in unsteady aerodynamics  Figure 21 is the uncertainty boundary of stochastic bifurcation diagram under two uncertainties bounds = 0.1 and = 0.2 for "Aerodynamic" model and "Structure" model. The figure shows that different parameters impact different position of bifurcation diagram; this is mainly because the uncertainty of parameters performs different sensitivities under different velocity. Figure 22 shows uncertainty boundary under two uncertainties bounds 1 = 0.1 and 1 = 0.2 for "Aerodynamic + Structure" model. Due to the superposition effect for aerodynamic uncertainty and structure uncertainty, the boundary can cover most of test data under 0.1 uncertainty bound; the boundary can cover all test data under 0.2 uncertainty. For getting an accuracy error boundary, some more error bound combinations were calculated; the result of one combination is shown in Figure 23. In this combination, uncertainty bound of aerodynamic parameters was set as 1 = 0.15, and uncertainty bound of structure parameters was set as 1 = 0.1. From the figure, it can be seen that the resulting boundary of this combination can totally cover the test data, which means that empirical aerodynamic parameters in the LB model are not so accurate; it may have 15% variation of its nominal value; structure model may have 10% variation of its nominal value.
International Journal of Aerospace Engineering

Conclusions
A deterministic aerodynamic model with modified LB model which is more accurate to low Mach number situation was built and verified by using test data. With the deterministic aerodynamic model, an aeroelastic model for airfoil with a third-order stiffness in both pitch and plunge DOF was built. Considering the stochastic parameters uncertainty for both aerodynamics and structure, the nonintrusive polynomial chaos expansion model to unsteady aerodynamics was utilized for the uncertain nonlinear aeroelastic analysis. From current work, the following can still be seen: (1) The modified LB model is more suitable for aerodynamics calculation at low Mach numbers. The nonintrusive polynomial chaos method can be applied to  Figure 20: Sensitivity of parameters in "Aerodynamic + Structure" model. aeroelastic limit cycle oscillations and dynamic stall flutter analysis with uncertainty. When the number of uncertainty parameters is small (less than 10), the PCE method has a higher calculation efficiency than MCS method. When a large number of uncertainty parameters is needed to be investigated, they can be put into groups separately.
(2) The uncertainty parameters have significant impact on normal force coefficient and moment coefficient in the aerodynamics model. The result boundary of normal distribution is tighter than the result boundary of uniform distribution. The four parameters and 1 , , and affect aerodynamic property most and were chosen to calculate the uncertainty of bifurcation.
(3) Aerodynamic uncertainty and structure uncertainty will impact different position of bifurcation diagram which is due to the fact that uncertainty parameters have different sensitivities under different velocity. Bifurcation is related to flow separation of aerodynamics, because the uncertainty of affects the bifurcation velocity most obviously in "Aerodynamic" model. Structure uncertainty also can impact the bifurcation position. Under same uncertainty bound, boundary variation for "Structure" model is larger than "Aerodynamic" model. (4) Taking both aerodynamic and structure uncertainty into consideration, when the uncertainty bound of aerodynamic parameters is 15% of its nominal value, at the same time, uncertainty bound of structure parameters is 10%; the range of uncertain LCO amplitude at all velocities can cover the experimental one. It indicates that the empirical aerodynamic parameters in the LB model are not so accurate; it may have 15% variation of its nominal value; structure model may have 10% variation of its nominal value.