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

Authors: Linpeng Wang; Yuting Dai; Chao Yang
Journal: International Journal of Aerospace Engineering
 (2017)

---

## Abstract

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.


---

## Body

## 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 [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]. 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 valueCN′′, which is the substitute value of CN′, is introduced.(1)CN′′=CN′-DbnDbn=Dbn-1exp⁡-ΔSTb+CNn′-CNn-1′exp⁡-ΔS2Tb,where Dbn 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ΔCNv is given as follows:(2)ΔCNv=B1f′′-fVx,where B1 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; Vx represents the shape function of normal force due to vortex, which is given by(3)Vx=sin3/2⁡πτ2Tv0<τ<Tvcos2⁡πτ-TvTvlτ>Tv,where τ is nondimensional time; Tv is time constant of vortex traveling over chord; Tvl is vortex passage time constant.The additional pitching moment coefficientΔCmv is also the effect of vortex convection over the airfoil. It is as follows:(4)ΔCmv=B21-cos⁡πτvTvlΔCNv,where B2 is the coefficient dependent on airfoil. τ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α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:(5)αmin=αmin⁡0+Trq,where Tr 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:(6)ΔCnrv=B1f′′-fVxr,where Vxr represents the shape function of normal force due to vortex at return process, which is given by(7)Vxr=sin3/2⁡πτr2Tr0<τr<Trcos2⁡πτr-Tr2Trτr>Tr,where τ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 Figure1. Considering the unsteady aerodynamic force, the motion equation of the wing section is written as(8)QhQθ=mSθSθIθh¨θ¨+Ch00Cθh˙θ˙+Kh00Kθhθ,where θ is the airfoil pitch angle. h is the vertical displacement. m is the mass per unit length. Kh and Kθ are spring stiffness in bending and torsion. Ch and Cθ are damping in bending and torsion. Sθ is static mass moment. Iθ is the polar moment of inertia about 1/4 chord. Qh and Qθ 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:(9)QhQθ=-12ρV2CLcl12ρV2CMc2l+12ρV2CLclxh,where CL and CM are lift coefficient and moment coefficient, respectively, c is the chord length, l is the span length, xh is the distance between 1/4 chord and elastic axis, ρ is air density, and V is the free stream velocity.Figure 1
Sketch of airfoil aeroelastic system.The angle of attackα and pitch ratio q which are used in LB model can be given as follows:(10)αq=θ+h˙Vθ˙cV.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 [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: CNα,α1,S1,S2,k0,k1,k2,CM0,η,CN1,Tp,Tf,Tvl,andTv and the two parameters in modified low Mach number LB model: Tb and Tr. All the uncertainties are noted as ξ=ξ1,ξ2,…,ξn; here n=16. Then responses in the time history are written as x(t,ξ), 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:(11)xt,ξ=∑j=1Px^jtψjξ,where ψjξ 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 ψj(ξ). (12)P=1+∑s=1d1s!∏r=0s-1n+r=n+d!n!d!.This is referred to as total-order expansion. In order to compute the stochastic output equation, it is required to solve a problem withP 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 μ and variance σ2. The statistics of the random solution are given by(13)μ=Ext,ξ=∫abxt,ξwξdξ=x^0tσ2=Varxt,ξ=∫abxt,ξ-x0t2wξdξ=∑i=1P-1xi2ψi2.With denoting the inner product(14)ψiψj≡∫wξψiξψjξdξ.It is possible to evaluate the deterministic functionsx^j(t) (j=0,1,…,P-1) by using the orthogonality of ψjξ.(15)x^jt=xt,ξψjξψj2ξ=∫abxt,ξψjξpξdξψj2ξ,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:(16)∫fξkρξkdξk≈Q1ikf≜∑jk=1n1ikfξk,jkikwk,jkik,where n1ik is the number of quadrature point ξk,jkik for one-dimensional quadrature rules. wk,jkik is the weight of quadrature point. Then, for multidimensional quadrature rules, the equation changes as follows:(17)Idf≜∫fξρξdξ≈Q1i1⊗⋯⊗Q1idf≜∑j1=1n1i1⋯∑jd=1n1idfξ1,j1i1,…,ξ1,jdjdw1,j1i1×⋯×wd,jdid.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(18)Le0=1,Le1=ξ,Lem+1=2m+1m+1ξLem-mm+1Lem-1.x ^ j ( t ) are deterministic functions to be evaluated and Lej(ξ) are stochastic part tabulated for uniform distribution. Then, (15) can be written as(19)x^jt=xt,ξLejξLej2=∫⋯∫fxξ1,…,ξdLejξ1,…,ξddξ1⋯dξdLej2.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.Table 1
Legendre polynomials quadrature point.
K Quadrature points w i 1 0 1 2 ± 0.5773502692 0.5 3 ± 0.7745966692 0.2777777778 0 0.4444444444Then, 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(20)CNα=CNα0+σ1ξ1α1=α10+σ2ξ2⋮Tr=Tr0+σnξn,where ξ1,ξ2,…,ξn are the random variables satisfying a uniform distribution and σ1,σ2,…,σn are their corresponding uncertainty bounds, respectively. Consequently, the response x in (11) can be rewritten as(21)x=∑i=0Pxiψξ1,ξ2,…,ξn.In (21), x0 represents the mean value of the response and xi is the corresponding fluctuation of the random components. Using the first-order representation of Legendre polynomials, the response can be written as (22)x=x0+ξ1x1+ξ2x2+⋯+ξ16x16.In the current work, normal aerodynamics coefficientCN and moment coefficient CM for modified low Mach number LB model and response of aeroelastic system are pitch angle of airfoil θ and its first-order derivative θ˙. They can be written as follows:(23)CN=x0,1+ξ1x1,1+ξ2x2,1+⋯+ξ16x16,1CM=x0,2+ξ1x1,2+ξ2x2,2+⋯+ξ16x16,2θ=x0,3+ξ1x1,3+ξ2x2,3+⋯+ξ16x16,3θ˙=x0,4+ξ1x1,4+ξ2x2,4+⋯+ξ16x16,4.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 (1/2π)e-x2/2 in the range of [-∞,+∞].The Legendre polynomial basis is written as(24)He0=1,He1=ξi,He2ξi1,ξi2=ξi12-1i1=i2ξi1ξi2i1≠i2.Then (15) becomes as follows:(25)x^jt=xt,ξHejξHej2=∫⋯∫fxξ1,…,ξdHejξ1,…,ξddξ1⋯dξdHej2.The information of quadrature point and its weight for Gauss–Hermite which is used for normal distribution is shown in Table2.Table 2
Hermit polynomials quadrature point.
K Quadrature points w i 1 0 1 2 ± 0.7071067812 0.5 3 ± 1.2247448714 0.1666666666667 0 0.6666666666667


## 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 valueCN′′, which is the substitute value of CN′, is introduced.(1)CN′′=CN′-DbnDbn=Dbn-1exp⁡-ΔSTb+CNn′-CNn-1′exp⁡-ΔS2Tb,where Dbn 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ΔCNv is given as follows:(2)ΔCNv=B1f′′-fVx,where B1 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; Vx represents the shape function of normal force due to vortex, which is given by(3)Vx=sin3/2⁡πτ2Tv0<τ<Tvcos2⁡πτ-TvTvlτ>Tv,where τ is nondimensional time; Tv is time constant of vortex traveling over chord; Tvl is vortex passage time constant.The additional pitching moment coefficientΔCmv is also the effect of vortex convection over the airfoil. It is as follows:(4)ΔCmv=B21-cos⁡πτvTvlΔCNv,where B2 is the coefficient dependent on airfoil. τ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α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:(5)αmin=αmin⁡0+Trq,where Tr 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:(6)ΔCnrv=B1f′′-fVxr,where Vxr represents the shape function of normal force due to vortex at return process, which is given by(7)Vxr=sin3/2⁡πτr2Tr0<τr<Trcos2⁡πτr-Tr2Trτr>Tr,where τr is a nondimensional time variable at return process.


## 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ΔCNv is given as follows:(2)ΔCNv=B1f′′-fVx,where B1 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; Vx represents the shape function of normal force due to vortex, which is given by(3)Vx=sin3/2⁡πτ2Tv0<τ<Tvcos2⁡πτ-TvTvlτ>Tv,where τ is nondimensional time; Tv is time constant of vortex traveling over chord; Tvl is vortex passage time constant.The additional pitching moment coefficientΔCmv is also the effect of vortex convection over the airfoil. It is as follows:(4)ΔCmv=B21-cos⁡πτvTvlΔCNv,where B2 is the coefficient dependent on airfoil. τ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α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:(5)αmin=αmin⁡0+Trq,where Tr 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:(6)ΔCnrv=B1f′′-fVxr,where Vxr represents the shape function of normal force due to vortex at return process, which is given by(7)Vxr=sin3/2⁡πτr2Tr0<τr<Trcos2⁡πτr-Tr2Trτr>Tr,where τ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 Figure1. Considering the unsteady aerodynamic force, the motion equation of the wing section is written as(8)QhQθ=mSθSθIθh¨θ¨+Ch00Cθh˙θ˙+Kh00Kθhθ,where θ is the airfoil pitch angle. h is the vertical displacement. m is the mass per unit length. Kh and Kθ are spring stiffness in bending and torsion. Ch and Cθ are damping in bending and torsion. Sθ is static mass moment. Iθ is the polar moment of inertia about 1/4 chord. Qh and Qθ 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:(9)QhQθ=-12ρV2CLcl12ρV2CMc2l+12ρV2CLclxh,where CL and CM are lift coefficient and moment coefficient, respectively, c is the chord length, l is the span length, xh is the distance between 1/4 chord and elastic axis, ρ is air density, and V is the free stream velocity.Figure 1
Sketch of airfoil aeroelastic system.The angle of attackα and pitch ratio q which are used in LB model can be given as follows:(10)αq=θ+h˙Vθ˙cV.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 [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: CNα,α1,S1,S2,k0,k1,k2,CM0,η,CN1,Tp,Tf,Tvl,andTv and the two parameters in modified low Mach number LB model: Tb and Tr. All the uncertainties are noted as ξ=ξ1,ξ2,…,ξn; here n=16. Then responses in the time history are written as x(t,ξ), 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:(11)xt,ξ=∑j=1Px^jtψjξ,where ψjξ 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 ψj(ξ). (12)P=1+∑s=1d1s!∏r=0s-1n+r=n+d!n!d!.This is referred to as total-order expansion. In order to compute the stochastic output equation, it is required to solve a problem withP 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 μ and variance σ2. The statistics of the random solution are given by(13)μ=Ext,ξ=∫abxt,ξwξdξ=x^0tσ2=Varxt,ξ=∫abxt,ξ-x0t2wξdξ=∑i=1P-1xi2ψi2.With denoting the inner product(14)ψiψj≡∫wξψiξψjξdξ.It is possible to evaluate the deterministic functionsx^j(t) (j=0,1,…,P-1) by using the orthogonality of ψjξ.(15)x^jt=xt,ξψjξψj2ξ=∫abxt,ξψjξpξdξψj2ξ,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:(16)∫fξkρξkdξk≈Q1ikf≜∑jk=1n1ikfξk,jkikwk,jkik,where n1ik is the number of quadrature point ξk,jkik for one-dimensional quadrature rules. wk,jkik is the weight of quadrature point. Then, for multidimensional quadrature rules, the equation changes as follows:(17)Idf≜∫fξρξdξ≈Q1i1⊗⋯⊗Q1idf≜∑j1=1n1i1⋯∑jd=1n1idfξ1,j1i1,…,ξ1,jdjdw1,j1i1×⋯×wd,jdid.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(18)Le0=1,Le1=ξ,Lem+1=2m+1m+1ξLem-mm+1Lem-1.x ^ j ( t ) are deterministic functions to be evaluated and Lej(ξ) are stochastic part tabulated for uniform distribution. Then, (15) can be written as(19)x^jt=xt,ξLejξLej2=∫⋯∫fxξ1,…,ξdLejξ1,…,ξddξ1⋯dξdLej2.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.Table 1
Legendre polynomials quadrature point.
K Quadrature points w i 1 0 1 2 ± 0.5773502692 0.5 3 ± 0.7745966692 0.2777777778 0 0.4444444444Then, 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(20)CNα=CNα0+σ1ξ1α1=α10+σ2ξ2⋮Tr=Tr0+σnξn,where ξ1,ξ2,…,ξn are the random variables satisfying a uniform distribution and σ1,σ2,…,σn are their corresponding uncertainty bounds, respectively. Consequently, the response x in (11) can be rewritten as(21)x=∑i=0Pxiψξ1,ξ2,…,ξn.In (21), x0 represents the mean value of the response and xi is the corresponding fluctuation of the random components. Using the first-order representation of Legendre polynomials, the response can be written as (22)x=x0+ξ1x1+ξ2x2+⋯+ξ16x16.In the current work, normal aerodynamics coefficientCN and moment coefficient CM for modified low Mach number LB model and response of aeroelastic system are pitch angle of airfoil θ and its first-order derivative θ˙. They can be written as follows:(23)CN=x0,1+ξ1x1,1+ξ2x2,1+⋯+ξ16x16,1CM=x0,2+ξ1x1,2+ξ2x2,2+⋯+ξ16x16,2θ=x0,3+ξ1x1,3+ξ2x2,3+⋯+ξ16x16,3θ˙=x0,4+ξ1x1,4+ξ2x2,4+⋯+ξ16x16,4.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 (1/2π)e-x2/2 in the range of [-∞,+∞].The Legendre polynomial basis is written as(24)He0=1,He1=ξi,He2ξi1,ξi2=ξi12-1i1=i2ξi1ξi2i1≠i2.Then (15) becomes as follows:(25)x^jt=xt,ξHejξHej2=∫⋯∫fxξ1,…,ξdHejξ1,…,ξddξ1⋯dξdHej2.The information of quadrature point and its weight for Gauss–Hermite which is used for normal distribution is shown in Table2.Table 2
Hermit polynomials quadrature point.
K Quadrature points w i 1 0 1 2 ± 0.7071067812 0.5 3 ± 1.2247448714 0.1666666666667 0 0.6666666666667


## 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 [27]. The freestream velocity is 0.1 Mach. The results at the condition of deep dynamic stall (angle of attack is α=15°+10°sin⁡ωt) are compared. The compared results are shown in Figures 2 and 3.Figure 2
Normal force coefficient atα=15°+10°sin⁡ωt.Figure 3
Moment coefficient atα=15°+10°sin⁡ωt.From Figures2 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 [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.Table 3
Important experimental parameters.
Parameter value C 0.3 m m 16.69 kg/m I θ 0.31 kgm2 k θ 13.1 Nm/rad l 0.9 m ρ 1.225 kg/m3 k h 30.5 N/mm Pitch axis 0.115 mTo 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(26)Kθθ+kθ3θ3,where kθ3 is the dimensionless parameter governing the third-order term. In DOF of plunge, similar expression of force can be given as (27)Khh+kh3h3.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.Figure 4
Pitch bifurcation diagram.Figure 5
Phase plane betweenθ and θ˙.From Figures4 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 θ and θ˙ 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.


## 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 [27]. The freestream velocity is 0.1 Mach. The results at the condition of deep dynamic stall (angle of attack is α=15°+10°sin⁡ωt) are compared. The compared results are shown in Figures 2 and 3.Figure 2
Normal force coefficient atα=15°+10°sin⁡ωt.Figure 3
Moment coefficient atα=15°+10°sin⁡ωt.From Figures2 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 [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.Table 3
Important experimental parameters.
Parameter value C 0.3 m m 16.69 kg/m I θ 0.31 kgm2 k θ 13.1 Nm/rad l 0.9 m ρ 1.225 kg/m3 k h 30.5 N/mm Pitch axis 0.115 mTo 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(26)Kθθ+kθ3θ3,where kθ3 is the dimensionless parameter governing the third-order term. In DOF of plunge, similar expression of force can be given as (27)Khh+kh3h3.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.Figure 4
Pitch bifurcation diagram.Figure 5
Phase plane betweenθ and θ˙.From Figures4 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 θ and θ˙ 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 isN(0,(0.0333)2) and N(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.Figure 6
Plot for different parameter distribution.
(a)
μ = 0, σ=0.0333 (normal), and σ1=0.1 (uniform) (b)
μ = 0, σ=0.0666 (normal), and σ1=0.2 (uniform)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. Figures7 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ξ 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 7
Mean of normal force coefficient in different corresponding uncertainty bounds for two kinds of distribution.Figure 8
Mean of moment coefficient in different corresponding uncertainty bounds for two kinds of distribution.Figure 9
Variance of normal force coefficientD(ξ) at upstroke and at return.
(a)
At upstroke (b)
At returnFigure10 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.Figure 10
Bound of MCS and nonintrusive PCE for two kinds of distribution.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 coefficientCN and moment coefficient CM for two kinds of distribution in (23) was gotten. Parameters were nondimensionalized by CN and CM separately. Figures 11 and 12 show uncertainty parameters in normal force coefficient CN for two kinds of distribution and Figure 13 shows uncertainty parameters in moment coefficient CM 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). CNα 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. Tp, Tvl, and Tvare more sensitive parameters for CN during upstroke than others, while, during return, Tf and Tr become more sensitive than other parameters. From Figure 13, it can be observed that parameters which are sensitive for CN are also sensitive for CM. However, the number of parameters which are sensitive for CM is more than those for CN, such as k0, k1, and k2. 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.Figure 11
Uncertainty parameters in normal force coefficientCN for normal distribution.
(a)
At upstroke (b)
At returnFigure 12
Uncertainty parameters in normal force coefficientCN for uniform distribution.
(a)
At upstroke (b)
At returnFigure 13
Uncertainty parameters in moment coefficientCM for uniform distribution.
(a)
At upstroke (b)
At returnFrom the parameters shown above, it can be seen that the parameters which are sensitive toCN also can effect CM obviously. For decreasing the calculating load, four parameters which have obvious effect on CN were filtered to use for investigating the stochastic uncertainty of LCO. The four parameters are CNα and α1, Tf, and Tr which will contribute most to normal for coefficient CN and moment coefficient CM. 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 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 stiffnessKθ and damping Cθ and spring stiffness Kh and damping Ch.Figures14, 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.Figure 14
Response of pitch angle atV=11 m/s.
(a)
Time history of pitch angle (b)
Phase plane betweenθ and θ˙Figure 15
Response of pitch angle atV=12.5 m/s.
(a)
Time history of pitch angle (b)
Phase plane betweenθ and θ˙Figure 16
Response of pitch angle atV=16 m/s.
(a)
Time history of pitch angle (b)
Phase plane betweenθ and θ˙Using the judgment rule defined above, bifurcation diagram with stochastic uncertainty was gotten. Figure17(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 V=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 17
Bifurcation diagram of pitch angle amplitude and sensitivity of parameters in “Aerodynamic.”
(a)
Stochastic bifurcation diagram (b)
Sensitivity of parametersFigure18(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 V=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 Kθ and Cθ were chosen as uncertainties parameters of structure to take part in building model “Aerodynamic + Structure” which contains both aerodynamic uncertainty and structure uncertainty.Figure 18
Bifurcation diagram of pitch angle amplitude and sensitivity of parameters in “Structure.”
(a)
Stochastic bifurcation diagram (b)
Sensitivity of parametersUsing 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. Figure19 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.Figure 19
Bifurcation diagram of pitch angle amplitude in “Aerodynamic + Structure.”Statistics information of bifurcation velocity for three models was shown in Table4. 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 σ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.Table 4
Statistics information of bifurcation velocity for three models.
Bifurcation PDF(0,1) (σ1=0.1) Bound variation(σ1=0.1) Bifurcation PDF(0,1) (σ1=0.2) Boundvariation(σ1=0.2) Bifurcation velocity (deterministic) Aerodynamic (11.5 m/s, 12.1 m/s) (−0.042, 0.008) (10.8 m/s, 12.1 m/s) (−0.100, 0.008) 12 m/s Structure (11.3 m/s, 12.3 m/s) (−0.058, 0.025) (10.8 m/s, 12.7 m/s) (−0.100, 0.058) Aerodynamic + Structure (11 m/s, 12.3 m/s) (−0.083, 0.025) (9.8 m/s, 12.7 m/s) (−0.183, 0.058)Figure20 is the sensitivities of parameters in “Aerodynamic + Structure”; it can be seen that structure uncertainty parameter Kθ is the most sensitive parameter for bifurcation; then aerodynamic uncertainty parameter Tf which is delay constant for the flow separation point due to unsteady effect is also a big sensitive parameter second only to Kθ. 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 20
Sensitivity of parameters in “Aerodynamic + Structure” model.Figure21 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 21
Uncertainty boundary in stochastic bifurcation diagram in “Aerodynamic”/“Structure.”
(a)
“Aerodynamic” (b)
“Structure”Figure22 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.Figure 22
Uncertainty boundary in stochastic bifurcation diagram in “Aerodynamic + Structure.”Figure 23
Uncertainty boundary for aerodynamic uncertainty bound 0.15 and structure uncertainty bound 0.1.


## 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 isN(0,(0.0333)2) and N(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.Figure 6
Plot for different parameter distribution.
(a)
μ = 0, σ=0.0333 (normal), and σ1=0.1 (uniform) (b)
μ = 0, σ=0.0666 (normal), and σ1=0.2 (uniform)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. Figures7 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ξ 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 7
Mean of normal force coefficient in different corresponding uncertainty bounds for two kinds of distribution.Figure 8
Mean of moment coefficient in different corresponding uncertainty bounds for two kinds of distribution.Figure 9
Variance of normal force coefficientD(ξ) at upstroke and at return.
(a)
At upstroke (b)
At returnFigure10 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.Figure 10
Bound of MCS and nonintrusive PCE for two kinds of distribution.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 coefficientCN and moment coefficient CM for two kinds of distribution in (23) was gotten. Parameters were nondimensionalized by CN and CM separately. Figures 11 and 12 show uncertainty parameters in normal force coefficient CN for two kinds of distribution and Figure 13 shows uncertainty parameters in moment coefficient CM 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). CNα 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. Tp, Tvl, and Tvare more sensitive parameters for CN during upstroke than others, while, during return, Tf and Tr become more sensitive than other parameters. From Figure 13, it can be observed that parameters which are sensitive for CN are also sensitive for CM. However, the number of parameters which are sensitive for CM is more than those for CN, such as k0, k1, and k2. 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.Figure 11
Uncertainty parameters in normal force coefficientCN for normal distribution.
(a)
At upstroke (b)
At returnFigure 12
Uncertainty parameters in normal force coefficientCN for uniform distribution.
(a)
At upstroke (b)
At returnFigure 13
Uncertainty parameters in moment coefficientCM for uniform distribution.
(a)
At upstroke (b)
At returnFrom the parameters shown above, it can be seen that the parameters which are sensitive toCN also can effect CM obviously. For decreasing the calculating load, four parameters which have obvious effect on CN were filtered to use for investigating the stochastic uncertainty of LCO. The four parameters are CNα and α1, Tf, and Tr which will contribute most to normal for coefficient CN and moment coefficient CM. 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 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 stiffnessKθ and damping Cθ and spring stiffness Kh and damping Ch.Figures14, 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.Figure 14
Response of pitch angle atV=11 m/s.
(a)
Time history of pitch angle (b)
Phase plane betweenθ and θ˙Figure 15
Response of pitch angle atV=12.5 m/s.
(a)
Time history of pitch angle (b)
Phase plane betweenθ and θ˙Figure 16
Response of pitch angle atV=16 m/s.
(a)
Time history of pitch angle (b)
Phase plane betweenθ and θ˙Using the judgment rule defined above, bifurcation diagram with stochastic uncertainty was gotten. Figure17(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 V=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 17
Bifurcation diagram of pitch angle amplitude and sensitivity of parameters in “Aerodynamic.”
(a)
Stochastic bifurcation diagram (b)
Sensitivity of parametersFigure18(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 V=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 Kθ and Cθ were chosen as uncertainties parameters of structure to take part in building model “Aerodynamic + Structure” which contains both aerodynamic uncertainty and structure uncertainty.Figure 18
Bifurcation diagram of pitch angle amplitude and sensitivity of parameters in “Structure.”
(a)
Stochastic bifurcation diagram (b)
Sensitivity of parametersUsing 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. Figure19 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.Figure 19
Bifurcation diagram of pitch angle amplitude in “Aerodynamic + Structure.”Statistics information of bifurcation velocity for three models was shown in Table4. 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 σ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.Table 4
Statistics information of bifurcation velocity for three models.
Bifurcation PDF(0,1) (σ1=0.1) Bound variation(σ1=0.1) Bifurcation PDF(0,1) (σ1=0.2) Boundvariation(σ1=0.2) Bifurcation velocity (deterministic) Aerodynamic (11.5 m/s, 12.1 m/s) (−0.042, 0.008) (10.8 m/s, 12.1 m/s) (−0.100, 0.008) 12 m/s Structure (11.3 m/s, 12.3 m/s) (−0.058, 0.025) (10.8 m/s, 12.7 m/s) (−0.100, 0.058) Aerodynamic + Structure (11 m/s, 12.3 m/s) (−0.083, 0.025) (9.8 m/s, 12.7 m/s) (−0.183, 0.058)Figure20 is the sensitivities of parameters in “Aerodynamic + Structure”; it can be seen that structure uncertainty parameter Kθ is the most sensitive parameter for bifurcation; then aerodynamic uncertainty parameter Tf which is delay constant for the flow separation point due to unsteady effect is also a big sensitive parameter second only to Kθ. 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 20
Sensitivity of parameters in “Aerodynamic + Structure” model.Figure21 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 21
Uncertainty boundary in stochastic bifurcation diagram in “Aerodynamic”/“Structure.”
(a)
“Aerodynamic” (b)
“Structure”Figure22 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.Figure 22
Uncertainty boundary in stochastic bifurcation diagram in “Aerodynamic + Structure.”Figure 23
Uncertainty boundary for aerodynamic uncertainty bound 0.15 and structure uncertainty bound 0.1.


## 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 parametersCNα and α1, Tf, and Tr 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 ofTf 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.



---

*Source: 2571253-2017-02-16.xml*