# Quantification of Limit Cycle Oscillations in Nonlinear Aeroelastic Systems with Stochastic Parameters

Authors: C. C. Cui; J. K. Liu; Y. M. Chen
Journal: International Journal of Aerospace Engineering
 (2016)

---

## Abstract

This paper presents an analytical feature for limit cycle oscillation (LCO) in the nonlinear aeroelastic system of an airfoil, with major emphasis on its applications in LCO quantification. The nonlinear stiffness is modeled as the product of thepth power of vibration displacement and the qth power of velocity, with its coefficient as a stochastic parameter. One interesting finding is that the LCO amplitude is directly proportional to the 1/(1-p-q)th power of the coefficient, whereas the frequency is independent of the coefficient. Based on this feature, the statistics and distribution functions of the LCO amplitude are obtained semianalytically, which are validated by Monte Carlo simulations. In addition, we discuss the possible influences of the nonlinear stiffness on flutter suppression of the airfoil subjected to Gaussian white noises. Surprisingly, increasing the nonlinear stiffness alone does not necessarily reduce the vibration amplitude as expected. Instead, it may sometimes induce disastrous subcritical LCOs with much higher vibration amplitudes.


---

## Body

## 1. Introduction

The nonlinear aeroelastic system of an airfoil is a typical self-excited system, which can exhibit lots of nonlinear dynamical behaviors such as bifurcations, limit cycle oscillations (LCO), and chaotic responses [1, 2]. Quantification of airfoil LCOs via analytical and/or semianalytical techniques has been an active area of research for many years. It has stimulated the interests and curiosities of many researchers [3–9]. Due to the design and manufacturing errors, uncertainties are usually inevitable in airfoil aeroelastic systems [10]. They usually happen to system parameters as stochastic [11] or uncertain-but-bounded variables [12]. As uncertainties are included, it will become much more cumbersome to investigate nonlinear aeroelastic responses such as LCOs, bifurcations, and chaos.Recent years have witnessed an increasing amount of research works focusing on the stochastic analysis of airfoil LCO. Expansion technique has been widely used to propose solution approaches on this topic. Mathematically, these methods aim at transforming a stochastic problem into a deterministic one, according to appropriate representations of the output variables or stochastic processes. Attar and Dowell [13] employed Monte Carlo simulation and a response surface method, respectively, to map the output characteristics of a delta wing into two input parameters. They considered two variables, that is, the thickness and modulus of elasticity. Millman et al. [14] also applied the response surface mapping to simulate the bifurcations of LCOs. Another technique, called stochastic Galerkin method, was proposed by Xiu et al. [15] to study LCOs arising in nonlinear systems with uncertainties. Pettit and Beran proposed methods for quantification of LCOs by using polynomial chaos [16] and Wiener-Haar expansion [17], respectively. Additionally, Millman et al. [18] applied Fourier chaos expansion to investigate the bifurcation behaviors of an airfoil. The Gegenbauer polynomial can also be applied in the orthogonal decomposition of stochastic nonlinear flutter system, as suggested by Wu et al. [19].Aside from the above-mentioned approaches based on expansion technique, there are some other approaches developed to solve nonlinear stochastic aeroelastic systems. The stochastic collocation method was applied by Deng et al. [20] to solve the nonlinear aeroelastic system of an airfoil with a control surface. In this method, deterministic solutions are obtained at a series of discrete values for stochastic parameters. Based on the attained solutions, the stochastic analysis can be carried out using an interpellation method. Chen et al. [21] applied the incremental harmonic method to solve a nonlinear aeroelastic system with uncertain-but-bounded parameters. In addition, the B-spline stochastic projection technique was utilized by Millman et al. to quantify uncertainties [22]. It was also extended, by Beran et al. [23], to the treatment of aerodynamic nonlinearities governed by discrete Euler equations. More recently, the normal form theory was adopted to investigate the flutter boundary of an airfoil [24].As mentioned above, a lot of numerical analysis has been made to better understand the effect of uncertain parameters on LCOs as well as other responses. Analytical or semianalytical prediction on aeroelastic responses is of fundamental significance to aeroelastics research; for example, Balakrishnan and Tuffaha established an explicit formula to determine aeroelastic divergence speed [25]. Nevertheless, it is rare to find more analytical and even semianalytical results. This study is motivated by an analytical dependence of the airfoil LCO on the nonlinear stiffness coefficient. More specifically, this parametric dependence regards the analytical relationship between the amplitude (frequency) and the coefficient of the nonlinear coefficient. It will be proved and extended, in Section 2, to general self-excited systems with multiple nonlinearities. Based on these features, we will present an algorithm to quantify the main statics of LCOs arising in the aeroelastic system with the nonlinear coefficient as a stochastic parameter. Furthermore, we will discuss the effects of the nonlinear coefficients on LCO suppression when the airfoil is subjected to Gaussian white noises.


## 2. Parametric Dependence of LCO

The nonlinear aeroelastic system [1, 2] of an airfoil is usually modeled as a set of ordinary differential equations. To elucidate the parametric dependence of the LCO, the system is rewritten as the following self-excited oscillator:(1)Mx¨+Cx˙+Kx+cFx,x˙=0in which the superscript denotes the differentiation with respect to time t, x is the unknown vector of dimension n, and M, C, and K are the mass, damping, and stiffness matrixes, respectively. The coefficient, c, is a nonzero constant. The nonlinear term, Fx,x˙, is given as a continuous function in one component of x:(2)Fx,x˙=0,…,0,xipx˙iq,0,…,0Twith p and q as given constants. The superscript “T” denotes the transpose and xi the ith component of vector x. The major result regarding the relationship between LCO and the nonlinear coefficient is given as follows:(I)
The LCO amplitude is directly proportional toc1/(1-p-q).(II)
The LCO frequency is independent ofc.Denote the LCO solution asx=y when c=1, with the angular frequency as ω. Under the transformation τ=ωt, (1) becomes(3)ω2Mx′′+ωCx′+Kx+cFx,ωx′=0,where the superscript denotes the differentiation with respect to τ. Considering that y is a solution of (3) as c=1, we have(4)ω2My′′+ωCy′+Ky+Fy,ωy′=0.For any nonzero constant λ, substituting x=λy into the left side of (3), we obtain(5)ω2My′′+ωCy′+Ky+λp+q-1cFy,ωy′.Comparing (4) and (5), we know that as long as λp+q-1c=1, x=λy is also a solution to (3); that is,(6)λ=c1/1-p-q.As the coefficient c varies, one of the LCOs can be expressed as c1/(1-p-q)y, which implies that the amplitude is in direct proportion to c1/(1-p-q) and that the frequency is independent of c.We take the famous van der Pol equation as an illustration(7)x¨+x-εx˙+δx2x˙=0.As is well known, as both ε and δ are positive constants, this system has a stable LCO. In this example, we have p=2 and q=1. According to the above properties, the amplitude of the LCO is in inverse proportion to δ.We examine the above-mentioned properties with the help of some highly accurate solutions. Tables1 and 2 present the results obtained by Delamotte [26], Buonomo [27], and the homotopy analysis method (HAM) [28, 29], respectively. Here, Delamotte’s and Buonomo’s solutions are provided, respectively, to validate the HAM ones. With such excellent agreement between these solutions, it is reasonable to consider the HAM results to be precise and use it as a benchmark for comparison in examining the above-mentioned features. One can observe that the frequency remains the same as nonlinear coefficient δ varies. On the other hand, as δ increases (or decreases) by 100 times, the amplitudes decreases (or increases) 10 times. That means the amplitude is inversely proportional to 100=10.Table 1
The frequency and amplitude of the van der Pol equation withε=1, by Delamotte’s and 100th-order HAM solution, respectively.
Solution a ω δ = 1  Delamotte 2.0086198609 0.94295584744161 HAM 2.00861986087484 0.94295584744161  δ = 0.25  HAM 4.01723972174969 0.94295584744161  δ = 0.01  HAM 20.08619860874844 0.94295584744161Table 2
The frequency and amplitude of the van der Pol equation withε=0.3, by Buonomo’s and 100th-order HAM solution, respectively.
Solution a ω δ = 0.3  Buonomo Not given 0.994419844392168 HAM 2.000922385554212 0.994419844392168  δ = 1.2  HAM 1.000461192777106 0.994419844392168  δ = 30  HAM 0.2000922385554212 0.994419844392168


## 3. Quantification of LCOs in Nonlinear Aeroelastic Systems

The physical model presented in Figure1 is a two-dimensional airfoil, oscillating in pitch and plunge, which has been employed by many authors [5–9].Figure 1
Sketch of a two-dimensional airfoil.b is the length of mid-chord, ab is the distance between the elastic axis and the mid-chord, and xb is the distance between the mass center and the elastic axis.The pitch angle about the elastic axis is denoted byα, positive with the nose up; the plunge deflection is denoted by h-, positive in the downward direction. In terms of nondimensional time t=ωαt1 (t1 is the real time) and nondimensional plunge displacement ξ=h-/b, the coupled motions of the airfoil in incompressible unsteady flow can be described as [30, 31](8)C0ξ¨+C1α¨+C2ξ˙+C3α˙+C4ξ+C5α+C6w1+C7w2+C8w3+C9w4+C10ξ+ηξ3=0,d0ξ¨+d1α¨+d2ξ˙+d3α˙+d4ξ+d5α+d6w1+d7w2+d8w3+d9w4+d10ξ+δα3=0,where η and δ are the coefficients of the cubic stiffness; and wi’s are given as (9)w1=∫0te-ε1t-σασdσ,w2=∫0te-ε2t-σασdσ,w3=∫0te-ε1t-σξσdσ,w4=∫0te-ε2t-σξσdσ,in which ε1=0.335 and ε2=0.0455. The coefficients and more details about the model are stated in Li et al. [30] or Liu and Dowell [31]. As one of the nonlinear stiffness coefficients (η and δ) is given as 0, the system contains a single nonlinearity. The nonlinear coefficients are given as bounded random parameters(10)η=η0+η1v,δ=δ0+δ1v,where η0, δ0 and η1, δ1 are positive constants and η0>η1, δ0>δ1; v∈[-1,1] is given as a stochastic parameter with the probability function as f(v).Case 1 (pitch nonlinearityη=0, δ≠0).
The system parameters in (8) are given in Table 3, and the value of nondimensional wind speed U varies. Numerical solutions are obtained by the Runge-Kutta method. Using analytical techniques developed for linear flutter analysis, the critical flutter speed is found to be U=Uf=6.0385. As U increases beyond Uf, a stable LCO arises. As U/Uf increases further, a secondary Hopf bifurcation arises following a jump for LCO amplitude [31].Table 3
Parameter values in system (8).
Parameter c 0 c 1 c 2 c 3 c 4 c 5 c 6 Value 1.01 0.255 0.01 0.02 0.0022 0.0122 1.4 E - 4  Parameter d 0 d 1 d 2 d 3 d 4 d 5 d 6 Value 1.02 1.015 0 0.04 0 0 0  Parameter c 7 c 8 c 9 c 10 d 7 d 8 d 9 d 10 Value 0.0014 - 6.8 E - 6 - 6 E - 8 0.0625U-2 0 0 0 U - 2We now turn to calculate the statistics of the LCO amplitudes. Since we consider a cubic nonlinearity such thatp=3 and q=0, the extremes of the LCOs, denoted as H for ξ and A for α, respectively, can be expressed as(11)H=H1δ,A=A1δ,where H(1) and A(1) are the extremes as δ (or η) equals 1. According to (11), the statistics and distribution functions of A and H can be determined by the following procedures. Note that this special case was revealed by Chen and Liu using the homotopy analysis method [32].Mathematically, the mean and standard deviation ofA can be, respectively, given as(12)
E
A
=
∫
-
1
1
A
δ
f
v
d
v
=
A
1
∫
-
1
1
f
v
δ
0
+
δ
1
v
d
v
,
(13)
D
A
=
∫
-
1
1
A
2
δ
f
v
d
v
-
E
A
2
=
A
2
1
∫
-
1
1
f
v
δ
0
+
δ
1
v
d
v
-
E
A
2
.A bounded probability density function with monopeak and symmetrically distributed within[-1,1], named as arc-like probability density function, is defined as follows [10]:(14)fv=21-v2π,v≤10,v>1.Since the amplitude of the LCO can be semianalytically expressed as monotonous function of δ, the minimum (Amin) and maximum (Amax) of A can be easily determined as δ varies in a bounded region. Moreover, the probable density function of Aδ can be obtained according to the probability theory as(15)ga=fδaδ′a,Amin≤a≤Amax.And, the distribution function can then be integrated as(16)Ga=∫Aminagada,Amin≤a≤Amax.Table4 shows the comparisons between the solutions provided by the presented method (i.e., (12) and (13)) and Monte Carlo simulation (MCS) results. The relative discrepancies are less than 0.5% when 5000 samples are taken into account. The distribution functions of A and H are plotted in Figures 2 and 3, respectively. Excellent agreement can also be observed. Note that the distribution curves for A and H are similar to each other. That is because the functions of A and H can be expressed by the same functions in δ.Table 4
Comparisons between the mean and standard derivation by the presented method and by MCS.
EA DA EH DH U = 1.5 U f, δ0=100, δ1=10  MCS, 5000 samples 0.1230 3.0929 E - 3 0.3179 8.0184 E - 3 Presented method 0.1233 3.0927 E - 3 0.3195 8.0945 E - 3 Relative difference 0.2% 0.007% 0.5% 0.09%  U = 2 U f, δ0=100, δ1=50  MCS, 5000 samples 0.2003 0.0271 0.6342 0.0858 Presented method 0.2009 0.0272 0.6315 0.0856 Relative difference 0.3% 0.4% 0.4% 0.2%Figure 2
Distribution functions ofA obtained by the presented method and by MCS, respectively. Parameters are given as U=2Uf, δ0=100, and δ1=50.Figure 3
Distribution functions ofH obtained by the presented method and by MCS, respectively. Parameters are given as U/Uf=2, δ0=100, and δ1=50.AsU/Uf increases beyond 2.1, there are three extremes for pitch displacement, as revealed by Liu and Dowell [31]. Table 5 shows the results obtained by the presented method for the means and derivations. They are both in excellent agreement with the MCS results. The distribution functions for the extremes are plotted in Figure 4. It also shows that the results of the presented method agree well with the MCS ones.Table 5
The mean and standard derivation ofA of the three extremes obtained by the presented method and by MCS with 5000 samples, respectively. Parameters are given as U/Uf=2.5, δ0=100, and δ1=50.
EA DA MCS 0.2881 0.0390 Presented method 0.2880 0.0390  MCS 0.1839 0.0249 Presented method 0.1839 0.0249  MCS 0.0759 0.0103 Presented method 0.0760 0.0103Figure 4
The distribution functions for the three extremes ofA with U/Uf=2.5, δ0=100, and δ1=50.Case 2 (plunge nonlinearityη≠0, δ=0).
As a cubic plunge stiffness is considered, the bifurcation chart obtained by the incremental harmonic balance method is plotted in Figure5. Note that the critical flutter speed depends on the linear coefficients rather than nonlinear coefficients. Different from the first case, a subcritical Hopf bifurcation is detected at Uf. As U/Uf decreases from 1 to 0.68 or so, there exists an unstable LCO. The LCO gains its stability at U/Uf=0.68, and a stable LCO always arises again as U/Uf increases from 0.68.Figure 5
Bifurcation chart of the aeroelastic system withδ=0 and η=80.Table6 presents the comparisons of the results by the presented method and MCS ones. The distribution functions are also obtained for both the pitch and plunge amplitudes, respectively. Excellent agreement can be also observed in Figure 6.Table 6
Comparisons of the mean and standard derivation by the presented method with those by MCS.
EA DA EH DH U = 1 U f, δ0=100, δ1=80  MCS, 5000 samples 1.2605 0.3273 0.7330 0.1902 Presented method 1.2646 0.3295 0.7339 0.1912 Relative difference 0.3% 0.7% 0.1% 0.5%  U = 2 U f, δ0=100, δ1=80  MCS, 4000 samples 1.3349 0.3494 1.0236 0.2640 Presented method 1.3364 0.3482 1.0269 0.2676 Relative difference 0.1% 0.3% 0.3% 1.3%Figure 6
The distribution functions for the three extremes ofA and H with U/Uf=1, δ0=100, and δ1=80.


## 4. Influence of Nonlinear Stiffness to LCO Suppression

For the aeroelastic system with a cubic stiffness, the LCO amplitude is inversely proportional to the square root of the nonlinear coefficient. Accordingly, one would like to reduce the LCO amplitude by raising the nonlinear stiffness. As an illustrative example, Figure7 shows that the amplitude reduces to a large extent as the nonlinear coefficient increases from 8 to 100. The initial condition (IC) is defined as [α0,α˙0,ξ0,ξ˙0]. The single fixed point of system (8) is the origin, which is unstable as U is larger than Uf. In such a case, the motions converge to the stable LCO regardless of ICs. As expected, increasing the cubic stiffness can suppress the LCO for the case of U/Uf>1.Figure 7
Time histories for system (8) with IC=[1°,0,0,0]. Parameters are U/Uf=1.5 and η=8 (solid lines) or η=100 (dashed lines).WhenU/Uf<1, however, increasing the cubic stiffness does not necessarily lead to decreasing of the LCO amplitude as expected. According to Figure 5, two LCOs coexist when U/Uf is located between 0.68 and 1, with the lower one being unstable and the other stable. Note that the single fixed point is also stable. That implies the motions may converge to either the fixed point or the stable LCO. By choosing different ICs, the convergent region of the LCO can be determined according to time histories obtained by time marching integration. As Figure 8 shows, there is critical boundary in the α0-h(0) plane (when α˙0=ξ˙0=0). The vibration approaches the stable LCO when the IC is given within the shadow region; otherwise it diminishes to the fixed point.Figure 8
The convergent region for the stable LCO and fixed point of system (8) with U/Uf=0.7; shadow region denotes convergent points for the LCO while the blank for the fixed point.Generally, it is useful for vibration suppression to enlarge the convergent region of the fixed point because it implies that the motions being convergent as noises happen is more probable. As mentioned above, the smaller the nonlinear coefficient is, the larger the convergent region of the fixed point becomes. Figure9 shows the time histories of system (8) with the same value for U but with different values for the nonlinear coefficient. As η=8, the motion converges to the fixed point. If η=100 is chosen, the motion approaches a LCO finally. In this case, increasing the cubic stiffness results in larger amplitude vibration rather than vibration suppression.Figure 9
Time histories for system (8) with IC=[1°,0,0,0]. Parameters are U/Uf=0.8 and η=8 (solid lines) or η=100 (dashed lines).Consider Gaussian white noises in both the pitching and plunging DOFs, such as(17)C0ξ¨+C1α¨+C2ξ˙+C3α˙+C4ξ+C5α+C6w1+C7w2+C8w3+C9w4+C10ξ+ηξ3=0.1ε1t,d0ξ¨+d1α¨+d2ξ˙+d3α˙+d4ξ+d5α+d6w1+d7w2+d8w3+d9w4+d10ξ+δα3=0.1ε2tin which the perturbations εi(t)’s are given as white Gaussian noises with mean values as 1. The time histories, as plotted in Figure 10, show that the vibration amplitude for η=8 is smaller than that for η=100. Actually, the vibrations resulted from the white noises as η=8. The vibration is too small to cause the vibration to reach the convergent region of the stable LCO. For the case of η=100, on the other hand, the motions are in essence LCOs subjected to noises.Figure 10
Time histories for system (17) with IC=[1°,0,0,0]. Parameters are U/Uf=0.8 and η=8 (solid lines) or η=100 (dashed lines).


## 5. Conclusions

We have investigated the LCO of an airfoil aeroelastic system with a nonlinear stiffness. The nonlinearity is expressed as a product ofpth power of vibration displacement and qth power of the velocity. Interestingly, the LCO amplitude is directly proportional to 1/(1-p-q)th power of the nonlinear coefficient, while the frequency is independent of this coefficient.Special emphasis has been placed on the applications of the above-mentioned feature. An illustrative application is in LCO quantification for the nonlinear aeroelastic system with the nonlinear coefficient as a stochastic parameter. The statistics and distribution functions of the LCO amplitudes can be semianalytically determined. Excellent agreement between the obtained results and the MCS ones validates the feasibility and efficiency of the applications in LCO quantification.In addition, the influence of the nonlinear stiffness on LCO suppression has been discussed in detail. When a stable LCO exists and the fixed point is unstable, that is, for the case when the wind speed increases beyond the critical flutter speed, the LCO can be suppressed by means of increasing the cubic stiffness. Below the flutter speed, however, increasing the cubic stiffness cannot reduce (as expected) but possibly result in larger magnitude of vibration.As is well known, parameters in nonlinear dynamical systems are closely related to system stability, bifurcations, chaos, and so on. Therefore, parametric study is one of major tasks in nonlinear dynamics analysis. This study could also be considered to be a preliminary and useful attempt at the dependence of vibration on system parameters (especially nonlinear coefficients) of self-excited systems. The simple yet efficient applications imply that the feature revealed in this paper could be more applicable in more nonlinear systems.



---

*Source: 3049163-2016-12-20.xml*