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

1. 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 . There are other methods to calculate nonlinear aerodynamics such as Leishman-Beddoes (LB) . 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 . Dai and Yang Reviewed the methods applied to aeroelastic uncertainty quantification , 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  and probabilistic one . 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 . Some related work has been carried out, including uncertainty quantification in aerodynamics and stochastic basis construction, Badcock et al. [16-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.

2. Uncertain LCO Analysis by the Nonintrusive PCE Method

2.1. 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 [C".sub.N], which is the substitute value of [C'.sub.N], is introduced.

[mathematical expression not reproducible] (1)

where [mathematical expression not reproducible] is an intermediate variable, Tb is the time lag constant, and S is the distance traveled by airfoil in semichords.

2.1.1. 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 [DELTA][C.sup.v.sub.N] is given as follows:

[DELTA][C.sup.v.sub.N] = [B.sub.1] (f" - f)[V.sub.x], (2)

where [B.sub.1] is a coefficient related to airfoils; f is separation location in terms of chord; f" is delayed separation location of f in the original LB model; [V.sub.x] represents the shape function of normal force due to vortex, which is given by

[mathematical expression not reproducible] (3)

where [tau] is nondimensional time; [T.sub.v] is time constant of vortex traveling over chord; [T.sub.vl] is vortex passage time constant.

The additional pitching moment coefficient [DELTA][C.sup.v.sub.m] is also the effect of vortex convection over the airfoil. It is as follows:

[mathematical expression not reproducible] (4)

where [B.sub.2] is the coefficient dependent on airfoil. [[tau].sub.v] is nondimensional time during vortex passage.

2.1.2. Modification of Normal Force and Pitching Moment Coefficients during Return. During return, there is also an overshoot for normal force coefficient. A value of [[alpha].sub.min 0] was given. It is assumed that it is the start of the convection process of overshoot. A value of [[alpha].sub.min] was given to define the end of the convection process of overshoot. Their relationship is as follows:

[[alpha].sub.min] = [[alpha].sub.min 0] + [T.sub.r]q, (5)

where [T.sub.r] is delay constant for reattachment process and q is reduced pitch rate.

The overshoot for normal force coefficient at return process is given as follows:

[DELTA][C.sup.v.sub.nr] = [B.sub.1] (f" - f)[V.sub.xr], (6)

where [V.sub.xr] represents the shape function of normal force due to vortex at return process, which is given by

[mathematical expression not reproducible] (7)

where [[tau].sub.r] is a nondimensional time variable at return process.

2.2. 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

[mathematical expression not reproducible] (8)

where [theta] is the airfoil pitch angle. h is the vertical displacement. m is the mass per unit length. [K.sub.h] and [K.sub.[theta]] are spring stiffness in bending and torsion. [C.sub.h] and [C.sub.[theta]] are damping in bending and torsion. [S.sub.[theta]] is static mass moment. [I.sub.[theta]] is the polar moment of inertia about 1/4 chord. [Q.sub.h] and [Q.sub.[theta]] 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:

[mathematical expression not reproducible] (9)

where [C.sub.L] and [C.sub.M] are lift coefficient and moment coefficient, respectively, c is the chord length, I is the span length, [x.sub.h] is the distance between 1/4 chord and elastic axis, [rho] is air density, and V is the free stream velocity.

The angle of attack a and pitch ratio q which are used in LB model can be given as follows:

[mathematical expression not reproducible] (10)

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.

2.3. 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 . 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: [mathematical expression not reproducible], [[alpha].sub.1], [S.sub.1], [S.sub.2], [k.sub.0], [k.sub.1], [k.sub.2], [mathematical expression not reproducible], [eta], [mathematical expression not reproducible], [T.sub.p], [T.sub.f], [T.sub.vl], and [T.sub.v] and the two parameters in modified low Mach number LB model: [T.sub.b] and [T.sub.r]. All the uncertainties are noted as [xi] = [[[xi].sub.1], [[xi].sub.2], ..., [[xi].sub.n]]; here n = 16. Then responses in the time history are written as x(t, [xi]), which will be also a stochastic vector.

PCE estimates its coefficients using the approaches random sampling or tensor-product quadrature . 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:

[mathematical expression not reproducible] (11)

where [[psi].sub.j]([xi]) are multidimensional orthogonal polynomial basis which are stochastic part. The expansion is truncated to P term which is determined by the number of variables n and the order d of polynomial [[psi].sub.j]([xi]).

[mathematical expression not reproducible] (12)

This is referred to as total-order expansion. In order to compute the stochastic output equation, it is required to solve a problem with P unknowns and P equations. Therefore, it is required to run the model for P times in order to evaluate the deterministic functions which are used to evaluate the mean [mu] and variance [[sigma].sup.2]. The statistics of the random solution are given by

[mathematical expression not reproducible] (13)

With ( ) denoting the inner product

([[psi].sub.i][[psi].sub.j]) = [integral] w ([xi]) [[psi].sub.i] ([xi]) [[psi].sub.j] ([xi]) d [xi]. (14)

It is possible to evaluate the deterministic functions [x.sub.j](t) (j = 0, 1, ..., P - 1)by using the orthogonality of [[psi].sub.j]([xi]).

[mathematical expression not reproducible] (15)

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:

[mathematical expression not reproducible] (16)

where [mathematical expression not reproducible] is the number of quadrature point [mathematical expression not reproducible] for one-dimensional quadrature rules. [mathematical expression not reproducible] is the weight of quadrature point. Then, for multidimensional quadrature rules, the equation changes as follows:

[mathematical expression not reproducible] (17)

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

[mathematical expression not reproducible] (18)

[x.sub.j](t) are deterministic functions to be evaluated and [Le.sub.j] ([xi]) are stochastic part tabulated for uniform distribution. Then, (15) can be written as

[mathematical expression not reproducible] (19)

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. K is the number of quadrature points, which can be gotten by using K = d + 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

[mathematical expression not reproducible] (20)

where [[xi].sub.1], [[xi].sub.2], ..., [[xi].sub.n] are the random variables satisfying a uniform distribution and [[sigma].sub.1], [[sigma].sub.2], ..., [[sigma].sub.n] are their corresponding uncertainty bounds, respectively. Consequently, the response x in (11) can be rewritten as

[mathematical expression not reproducible] (21)

In (21), [x.sub.0] represents the mean value of the response and [x.sub.i] is the corresponding fluctuation of the random components. Using the first-order representation of Legendre polynomials, the response can be written as

x = [x.sub.0] + [[xi].sub.1][x.sub.1] + [[xi].sub.2][x.sub.2] + ... + [[xi].sub.16][x.sub.16]. (22)

In the current work, normal aerodynamics coefficient [C.sub.N] and moment coefficient [C.sub.M] for modified low Mach number LB model and response of aeroelastic system are pitch angle of airfoil [theta] and its first-order derivative [theta]. They can be written as follows:

[mathematical expression not reproducible] (23)

For calculating normal distribution, the Hermit polynomials are applied. In this polynomial series, the range of random variables is located in [-[infinity], +[infinity]]. The probability density function of each variable is [mathematical expression not reproducible] in the range of [-[infinity], +[infinity]].

The Legendre polynomial basis is written as

[mathematical expression not reproducible] (24)

Then (15) becomes as follows:

[mathematical expression not reproducible] (25)

The information of quadrature point and its weight for Gauss-Hermite which is used for normal distribution is shown in Table 2.

3. 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.

3.1. 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 . The freestream velocity is 0.1 Mach. The results at the condition of deep dynamic stall (angle of attack is [alpha] = 15[degrees] + 10[degrees] sin wt) 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.

3.2. Validation for LCO Analysis. Reference  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 a nonlinear incarnation of the standard symmetric airfoil with pitch and plunge. The original deterministic aeroelastic 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

[mathematical expression not reproducible] (26)

where [mathematical expression not reproducible] is the dimensionless parameter governing the third-order term. In DOF of plunge, similar expression of force can be given as

[mathematical expression not reproducible] (27)

The result of present work was compared with the result of experiment and the result of Song , 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 [theta] and [theta] 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.

4. 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.

4.1. 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 N(0, [(0.0333).sup.2]) and N(0, [(0.0666).sup.2]), the mean [mu] = 0, and standard deviation is [sigma] = 0.0333 and [sigma] = 0.0666. The corresponding uncertainty bounds of aerodynamic parameters for uniform distribution are set as [[sigma].sub.1] = 0.1 and [[sigma].sub.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 D([xi]) 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 ([sigma] = 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 [C.sub.N] and moment coefficient [C.sub.M] for two kinds of distribution in (23) was gotten. Parameters were nondimensionalized by [C.sub.N] and [C.sub.M] separately. Figures 11 and 12 show uncertainty parameters in normal force coefficient [C.sub.N] for two kinds of distribution and Figure 13 shows uncertainty parameters in moment coefficient [C.sub.M] 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). [mathematical expression not reproducible] and [[alpha].sub.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. [T.sub.p], [T.sub.vl], and [T.sub.v] are more sensitive parameters for [C.sub.N] during upstroke than others, while, during return, [T.sub.f] and [T.sub.r] become more sensitive than other parameters. From Figure 13, it can be observed that parameters which are sensitive for [C.sub.N] are also sensitive for [C.sub.M]. However, the number of parameters which are sensitive for CM is more than those for [C.sub.N], such as [k.sub.0], [k.sub.1], and [k.sub.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 [C.sub.N] also can effect [C.sub.M] obviously. For decreasing the calculating load, four parameters which have obvious effect on [C.sub.N] were filtered to use for investigating the stochastic uncertainty of LCO. The four parameters are [mathematical expression not reproducible] and [[alpha].sub.1], [T.sub.f], and [T.sub.r] which will contribute most to normal for coefficient [C.sub.N] and moment coefficient [C.sub.M]. 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 is similar. Based on the reasons above, when investigating uncertainty of LCO and bifurcation, uniform distribution of parameters was used in calculation.

4.2. 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 onlyaerodynamic 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 [K.sub.[theta]] and damping [C.sub.[theta]] and spring stiffness [K.sub.h] and damping [C.sub.h].

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 V = 16 m/s is larger than the one at V = 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 [sigma].sub.1] = 0.1 and [sigma].sub.1] = 0.2 in model "Aerodynamic" which are compared with the deterministic value and Figure 17(b) shows sensitivity of parameters for [sigma].sub.1] = 0.2 and probability density function (PDF) of bifurcation. In order to explain the problem clearly, the parameters were nondimensionalized by divided peak [theta] 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 [sigma].sub.1] = 0.2 is V - 10.8 m/s which is smaller than the mean value of 11.5 m/s with the corresponding uncertainty bound [sigma] = 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 [[sigma].sub.1] = 0.1 and [sigma].sub.1] = 0.2 in model "Structure" which are compared with the deterministic value and Figure 18(b) shows sensitivities of parameters for [sigma].sub.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 [sigma].sub.1] = 0.2 is V = 10.8 m/s which is smaller than the mean value of 11.3 m/s with the corresponding uncertainty bound [sigma] = 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 [K.sub.[theta]] and [C.sub.[theta]] 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 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 [[sigma].sub.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 [sigma].sub.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 [sigma].sub.1] = 0.1, the bound variation is 0.108; when [sigma].sub.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 [K.sub.[theta]] is the most sensitive parameter for bifurcation; then aerodynamic uncertainty parameter [T.sub.f] which is delay constant for the flow separation point due to unsteady effect is also a big sensitive parameter second only to [K.sub.[theta]]. The delay constant of the separation point in unsteady aerodynamics makes a great effect on the bifurcation point and it also affects the amplitude of the pitch angle.

Figure 21 is the uncertainty boundary of stochastic bifurcation diagram under two uncertainties bounds [sigma] = 0.l and [sigma] = 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 [sigma].sub.1] = 0.l and [sigma].sub.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 [sigma].sub.1] = 0.l5, and uncertainty bound of structure parameters was set as [sigma].sub.1] = 0.l. 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.

5. 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 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 [mathematical expression not reproducible] and [[alpha].sub.1], [T.sub.f], and [T.sub.r] 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 [T.sub.f] 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.

http://dx.doi.org/10.1155/2017/2571253

Competing Interests

The authors declare that they have no competing interests.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (nos. 11302011, 11402013) and the Research Fund for the Doctoral Program of Higher Education of China (no. 20131102120051).

References

 B. Timothy, "Non-intrusive uncertainty propagation with error bounds for conservation laws containing discontinuities," Uncertainty Quantification in Computational Fluid Dynamics, vol. 92, pp. 1-57, 2013.

 J. G. Leishman and T. S. Beddoes, "A semi-empirical model for dynamic stall," Journal of the American Helicopter Society, vol. 34, no. 3, pp. 3-17, 1989.

 D. Ragni, B. W. van Oudheusden, and F. Scarano, "Non-intrusive aerodynamic loads analysis of an aircraft propeller blade," Experiments in Fluids, vol. 51, no. 2, pp. 361-371, 2011.

 Y. Gu, X. Zhang, and Z. Yang, "Robust flutter analysis based on genetic algorithm," Science China Technological Sciences, vol. 55, no. 9, pp. 2474-2481, 2012.

 D. M. Pitt, D. P. Haudrich, M. J. Thomas et al., "Probabilistic aeroelastic analysis and its implications on flutter margin requirements," AIAA-2008-2198, 2008.

 Y. Dai and C. Yang, "Methods and advances in the study of aeroelasticity with uncertainties," Chinese Journal of Aeronautics, vol. 27, no. 3, pp. 461-474, 2014.

 Z. Wu, Y. Dai, C. Yang, and L. Chen, "Aeroelastic wind-tunnel test for aerodynamic uncertainty model validation," Journal of Aircraft, vol. 50, no. 1, pp. 47-55, 2013.

 Y. Dai and C. Yang, "Aeroservoelastic model modification and uncertainty quantification with Bayesian Posteriori Estimation," Transactions of the Japan Society for Aeronautical and Space Sciences, vol. 58, no. 4, pp. 237-246, 2015.

 Y. Dai, Z. Wu, and C. Yang, "Identification and robust limit-cycle- oscillation analysis of uncertain aeroelastic system," Science China Technological Sciences, vol. 54, no. 7, pp. 1841-1848, 2011.

 D. Borglund and U. Ringertz, "Solution of the flutter eigen-value problem with mixed structural/aerodynamic uncertainty," Journal of Aircraft, vol. 48, no. 1, pp. 343-348, 2011.

 C. L. Pettit, "Uncertainty quantification in aeroelasticity: recent results and research challenges," Journal of Aircraft, vol. 41, no. 5, pp. 1217-1229, 2004.

 H. N. Najm, "Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics," Annual Review of Fluid Mechanics, vol. 41, pp. 35-52, 2009.

 L. Mehrez, A. Doostan, D. Moens, and D. Vandepitte, "Stochastic identification of composite material properties from limited experimental databases, part ii: uncertainty modelling," Mechanical Systems and Signal Processing, vol. 27, pp. 484-498, 2012.

 M. Nikbay and P. Acar, "Reduced order modelling for static and dynamic aeroelastic predictions with multidisciplinary approach," CEAS Aeronautical Journal, vol. 6, no. 3, pp. 455-469, 2015.

 P. S. Beran, C. L. Pettit, and D. R. Millman, "Uncertainty quantification of limit-cycle oscillations," Journal of Computational Physics, vol. 217, no. 1, pp. 217-247, 2006.

 K. J. Badcock, S. Timme, S. Marques et al., "Transonic aeroelastic simulation for instability searches and uncertainty analysis," Progress in Aerospace Sciences, vol. 47, no. 5, pp. 392-423, 2011.

 L. Bruno, C. Canuto, and D. Fransos, "Stochastic aerodynamics and aeroelasticity of a flat plate via generalised Polynomial Chaos," Journal of Fluids and Structures, vol. 25, no. 7, pp. 1158-1176, 2009.

 S. Song, Z. Lu, W. Zhang, and L. Cui, "Uncertainty importance measure by fast fourier transform for wing transonic flutter," Journal of Aircraft, vol. 48, no. 2, pp. 449-455, 2011.

 D. K. Kishor, R. Ganguli, and S. Gopalakrishnan, "Uncertainty analysis of vibrational frequencies of an incompressible liquid in a rectangular tank with and without a baffle using polynomial chaos expansion," Acta Mechanica, vol. 220, no. 1, pp. 257-273, 2011.

 Z. Qiu, R. Huang, X. Wang, and W. Qi, "Structural reliability analysis and reliability-based design optimization: recent advances," Science China: Physics, Mechanics and Astronomy, vol. 56, no. 9, pp. 1611-1618, 2013.

 C. L. Pettit and P. S. Beran, "Spectral and multiresolution Wiener expansions of oscillatory stochastic processes," Journal of Sound and Vibration, vol. 294, no. 4, pp. 752-779, 2006.

 C. L. Pettit and P. S. Beran, "Polynomial chaos expansion applied to airfoil limit cycle oscillations," in Proceedings of the 45thA- IAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, AIAA Paper 2004-1691, Palm Springs, Calif, USA, 2004.

 Y. Dai and C. Yang, "Smolyak-grid-based flutter analysis with the stochastic aerodynamic uncertainty," Discrete Dynamics in Nature and Society, vol. 2014, Article ID 174927, 8 pages, 2014.

 D. Xiu, D. Lucor, C.-H. Su, and G. E. Karniadakis, "Stochastic modeling of flow-structure interactions using generalized polynomial chaos," Journal of Fluids Engineering, vol. 124, no. 1, pp. 51-59, 2002.

 G. Remeo, G. Frulla, E. Cestino, and G. Corsino, "HELIPLAT: design, aerodynamic, structural design of UAV scaled prototype," in Proceedings of the Conference of the 17th Italian Association of Aeronautics and Astronautics National Congress (AIDAA '03), vol. 1, pp. 569-578, Roma, Italy, September 2003.

 G. Frulla and E. Cestino, "Design, manufacturing and testing of a HALE-UAV structural demonstrator," Composite Structures, vol. 83, no. 2, pp. 143-153, 2008.

 W. Sheng, R. A. M. Galbraith, and F. N. Coton, "A modified dynamic stall model for low mach numbers," Journal of Solar Energy Engineering, Transactions of the ASME, vol. 130, no. 3, Article ID 031013, 10 pages, 2008.

 J. G. Leishman and G. L. Crouse, "State-space model for unsteady airfoil behavior and dynamic stall," in Proceedings of the 30th Structures, Structural Dynamics and Materials Conference, Structures, Structural Dynamics, and Materials and Co-Located Conferences, AIAA 89-1319-CP, Mobile, Ala, USA, April 1989.

 D. Xiu and G. E. Karniadakis, "Modeling uncertainty in flow simulations via generalized polynomial chaos," Journal of Computational Physics, vol. 187, no. 1, pp. 137-167, 2003.

 G. Dimitriadis and J. Li, "Bifurcation behavior of airfoil undergoing stall flutter oscillations in low-speed wind tunnel," AIAA Journal, vol. 47, no. 11, pp. 2577-2596, 2009.

 S. Shao, Q. Zhu, C. Zhang, and X. Ni, "Airfoii aeroelastic flutter analysis based on modified Leishman-Beddoes model at low Mach number," Chinese Journal of Aeronautics, vol. 24, no. 5, pp. 550-557, 2011.

Linpeng Wang, Yuting Dai, and Chao Yang

School of Aeronautic Science and Engineering, Beihang University, Beijing 100191, China

Correspondence should be addressed to Yuting Dai; yutingdai@buaa.edu.cn

Received 18 October 2016; Revised 14 December 2016; Accepted 19 December 2016; Published 16 February 2017

Caption: Figure 1: Sketch of airfoil aeroelastic system.

Caption: Figure 2: Normal force coefficient at [alpha] = 15[degrees] + 10[degrees] sin [omega]t.

Caption: Figure 3: Moment coefficient at [alpha] = 15[degrees] + 10[degrees] sin [omega]t.

Caption: Figure 4: Pitch bifurcation diagram.

Caption: Figure 5: Phase plane between [theta] and [theta].

Caption: Figure 6: Plot for different parameter distribution.

Caption: Figure 7: Mean of normal force coefficient in different corresponding uncertainty bounds for two kinds of distribution.

Caption: Figure 8: Mean of moment coefficient in different corresponding uncertainty bounds for two kinds of distribution.

Caption: Figure 9: Variance of normal force coefficient D([xi]) at upstroke and at return.

Caption: Figure 10: Bound of MCS and nonintrusive PCE for two kinds of distribution.

Caption: Figure 11: Uncertainty parameters in normal force coefficient [C.sub.N] for normal distribution.

Caption: Figure 12: Uncertainty parameters in normal force coefficient [C.sub.N] for uniform distribution.

Caption: Figure 13: Uncertainty parameters in moment coefficient [C.sub.M] for uniform distribution.

Caption: Figure 14: Response of pitch angle at V = 11 m/s.

Caption: Figure 15: Response of pitch angle at V = 12.5 m/s.

Caption: Figure 16: Response of pitch angle at V = 16 m/s.

Caption: Figure 17: Bifurcation diagram of pitch angle amplitude and sensitivity of parameters in "Aerodynamic."

Caption: Figure 18: Bifurcation diagram of pitch angle amplitude and sensitivity of parameters in "Structure."

Caption: Figure 19: Bifurcation diagram of pitch angle amplitude in "Aerodynamic + Structure."

Caption: Figure 20: Sensitivity of parameters in "Aerodynamic + Structure" model.

Caption: Figure 21: Uncertainty boundary in stochastic bifurcation diagram in "Aerodynamic"/"Structure."

Caption: Figure 22: Uncertainty boundary in stochastic bifurcation diagram in "Aerodynamic + Structure."

Caption: Figure 23: Uncertainty boundary for aerodynamic uncertainty bound 0.15 and structure uncertainty bound 0.1.
```Table 1: Legendre polynomials quadrature point.

1              0                    1
2    [+ or -] 0.5773502692         0.5
3     [+ or -] 07745966692    0.2777777778
0              0.4444444444

Table 2: Hermit polynomials quadrature point.

1              0                     1
2    [+ or -] 0.7071067812          0.5
3     [+ or -] 12247448714    0.1666666666667
0              0.6666666666667

Table 3: Important experimental parameters.

Parameter value

C                        0.3 m
m                      16.69 kg/m
[I.sub.[theta]]     0.31 kg[m.sup.2]
l                        0.9 m
[rho]              1.225 kg/[m.sup.3]
[k.sub.h]              30.5 N/mm
Pitch axis              0.115 m

Table 4: Statistics information of bifurcation velocity for three
models.

Bifurcation PDF
(0,1)                 Bound variation
([[sigma].sub.1] = -0.1)   ([[sigma].sub.1] = 0.1)

Aerodynamic     (11.5 m/s, 12.1 m/s)         (-0.042, 0.008)
Structure       (11.3 m/s, 12.3 m/s)         (-0.058, 0.025)
Aerodynamic      (11 m/s, 12.3 m/s)          (-0.083, 0.025)
+ Structure

Bifurcation PDF
(0,1)                Bound variation
([[sigma].sub.1] = 0.2)   ([[sigma].sub.1] = 0.2)

Aerodynamic    (10.8 m/s, 12.1 m/s)         (-0.100, 0.008)
Structure      (10.8 m/s, 12.7 m/s)         (-0.100, 0.058)
Aerodynamic     (9.8 m/s, 12.7 m/s)         (-0.183, 0.058)
+ Structure

Bifurcation
velocity
(deterministic)

Aerodynamic
Structure         12 m/s
Aerodynamic
+ Structure
```
Title Annotation: Printer friendly Cite/link Email Feedback Research Article; polynomial chaos expansion Wang, Linpeng; Dai, Yuting; Yang, Chao International Journal of Aerospace Engineering Report 1USA Jan 1, 2017 7237 Experimental Design Validation of Tilting Calibration Mechanism by Using Shape Memory Alloy Spring Actuator. Efficient Energy Flight Path Planning Algorithm Using 3-D Visibility Roadmap for Small Unmanned Aerial Vehicle. Aerodynamics Aeroelasticity Monte Carlo method Monte Carlo methods