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

This paper presents an analytical feature for limit cycle oscillation (LCO) in the nonlinear aeroelastic system of an airfoil, withmajor emphasis on its applications in LCO quantification.The nonlinear stiffness is modeled as the product of the pth 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.


Introduction
The nonlinear aeroelastic system of an airfoil is a typical selfexcited 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][4][5][6][7][8][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-butbounded 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-butbounded 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.

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: in which the superscript denotes the differentiation with respect to time , x is the unknown vector of dimension , and M, C, and K are the mass, damping, and stiffness matrixes, respectively.The coefficient, , is a nonzero constant.The nonlinear term, F(x, ẋ ), is given as a continuous function in one component of x: with  and  as given constants.The superscript "T" denotes the transpose and x  the th 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 to  1/(1−−) .
(II) The LCO frequency is independent of .Denote the LCO solution as x = y when  = 1, with the angular frequency as .Under the transformation  = , (1) becomes where the superscript denotes the differentiation with respect to .Considering that y is a solution of (3) as  = 1, we have For any nonzero constant , substituting x = y into the left side of (3), we obtain Comparing ( 4) and ( 5), we know that as long as  +−1  = 1, x = y is also a solution to (3); that is, As the coefficient  varies, one of the LCOs can be expressed as  1/(1−−) y, which implies that the amplitude is in direct proportion to  1/(1−−) and that the frequency is independent of .
We take the famous van der Pol equation as an illustration As is well known, as both  and  are positive constants, this system has a stable LCO.In this example, we have  = 2 and  = 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.Tables 1 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 abovementioned 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.

Quantification of LCOs in Nonlinear Aeroelastic Systems
The physical model presented in Figure 1 is a two-dimensional airfoil, oscillating in pitch and plunge, which has been employed by many authors [5][6][7][8][9].
The pitch angle about the elastic axis is denoted by , positive with the nose up; the plunge deflection is denoted by ℎ, positive in the downward direction.In terms of nondimensional time  =    1 ( 1 is the real time) and nondimensional plunge displacement  = ℎ/, the coupled motions of the airfoil in incompressible unsteady flow can be described as [30,31] where  and  are the coefficients of the cubic stiffness; and   's are given as (9) 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 where  0 ,  0 and  1 ,  1 are positive constants and is given as a stochastic parameter with the probability function as (V).
Case 1 (pitch nonlinearity  = 0,  ̸ = 0).The system parameters in (8) are given in Table 3, and the value of nondimensional wind speed  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  =   = 6.0385.As  increases beyond   , a stable LCO arises.As /  increases further, a secondary Hopf bifurcation arises following a jump for LCO amplitude [31].
We now turn to calculate the statistics of the LCO amplitudes.Since we consider a cubic nonlinearity such that  = 3 and  = 0, the extremes of the LCOs, denoted as  for  and  for , respectively, can be expressed as where (1) and (1) are the extremes as  (or ) equals 1.
According to (11), the statistics and distribution functions of  and  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 of  can be, respectively, given as  A bounded probability density function with monopeak and symmetrically distributed within [−1, 1], named as arclike probability density function, is defined as follows [10]: Since the amplitude of the LCO can be semianalytically expressed as monotonous function of , the minimum ( min ) and maximum ( max ) of  can be easily determined as  varies in a bounded region.Moreover, the probable density function of () can be obtained according to the probability theory as And, the distribution function can then be integrated as Table 4 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  and  are plotted in Figures 2 and 3, respectively.Excellent agreement can also be observed.Note that the distribution curves for  and  are similar to each other.That is because the functions of  and  can be expressed by the same functions in .
As /  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.
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 Figure 5.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   .As /  decreases from 1 to 0.68 or so, there exists an unstable LCO.The LCO gains its stability at /  = 0.68, and a stable LCO always arises again as /  increases from 0.68.are also obtained for both the pitch and plunge amplitudes, respectively.Excellent agreement can be also observed in Figure 6.

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 A (degree) Presented method MCS, 5000 samples  an illustrative example, Figure 7 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  is larger than   .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 /  > 1.  /  < 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 /  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) − ℎ(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.
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.Figure 9 shows the time histories of system (8) with the same value for  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.
Consider Gaussian white noises in both the pitching and plunging DOFs, such as in which the perturbations   ()'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.

Conclusions
We have investigated the LCO of an airfoil aeroelastic system with a nonlinear stiffness.The nonlinearity is expressed as a product of th power of vibration displacement and th power of the velocity.Interestingly, the LCO amplitude is directly proportional to 1/(1 −  − )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.

Figure 1 :
Figure1: Sketch of a two-dimensional airfoil. is the length of midchord,  is the distance between the elastic axis and the mid-chord, and  is the distance between the mass center and the elastic axis.

Figure 8 :Figure 9 :
Figure8: The convergent region for the stable LCO and fixed point of system(8) with /  = 0.7; shadow region denotes convergent points for the LCO while the blank for the fixed point.

Table 1 :
The frequency and amplitude of the van der Pol equation with  = 1, by Delamotte's and 100th-order HAM solution, respectively.

Table 2 :
The frequency and amplitude of the van der Pol equation with  = 0.3, by Buonomo's and 100th-order HAM solution, respectively.

Table 4 :
Comparisons between the mean and standard derivation by the presented method and by MCS.

Table 5 :
The mean and standard derivation of A of the three extremes obtained by the presented method and by MCS with 5000 samples, respectively.Parameters are given as /  = 2.5,  0 = 100, and  1 = 50.

Table 6
presents the comparisons of the results by the presented method and MCS ones.The distribution functions Figure 2: Distribution functions of  obtained by the presented method and by MCS, respectively.Parameters are given as  = 2  ,  0 = 100, and  1 = 50.

Table 6 :
Comparisons of the mean and standard derivation by the presented method with those by MCS.