Limit cycle oscillation behavior of transonic control surface buzz considering free-play nonlinearity

Journal of Fluids and Structures 61 (2016) 431–449

Contents lists available at ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

Limit cycle oscillation behavior of transonic control surface buzz considering free-play nonlinearity Shun He, Zhichun Yang n, Yingsong Gu School of Aeronautics, Northwestern Polytechnical University, Xi'an 710072, People's Republic of China

a r t i c l e i n f o

abstract

Article history: Received 8 April 2015 Accepted 20 November 2015 Available online 21 January 2016

The limit cycle oscillation (LCO) behaviors of control surface buzz in transonic flow are studied. Euler equations are employed to obtain the unsteady aerodynamic forces for Type B and Type C buzz analyses, and an all-movable control surface model, a wing/control surface model and a three-dimensional wing with a full-span control surface are adopted in the study. Aerodynamic and structural describing functions are used to deal with aerodynamic and structural nonlinearities, respectively. Then the buzz speed and buzz frequency are obtained by V-g method. The LCO behavior of the transonic control surface buzz system with linear structure exhibits subcritical or supercritical bifurcation at different Mach numbers. For nonlinear structural model with a free-play nonlinearity in the control surface deflection stiffness, the double LCO phenomenon is observed in certain range of flutter speed. The free-play nonlinearity changes the stability of LCOs at small amplitudes and turns the unstable LCO into a stable one. The LCO behavior is dominated by the aerodynamic nonlinearity for the case with large control surface oscillation amplitude but by the structural nonlinearity for the case with small amplitude. Good agreements between LCO behaviors obtained by the present method and available experimental data show that our study may help to explain the experimental observation in wind tunnel tests and to understand the physical mechanism of transonic control surface buzz. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Control surface buzz Transonic flow Aerodynamic nonlinearity Free-play Single degree of freedom flutter Limit cycle oscillation

1. Introduction Buzz is an aeroelastic phenomenon that occurs on the aircraft control surface in transonic flow, which is treated as a kind of the single degree of freedom flutter. It was firstly observed in aircraft flight around 1945 according to Lambourne (1964). During the flight test of a P-80 jet fighter aircraft (Lambourne, 1964), a low amplitude control surface buzz was encountered, and in the following flight tests of the P-80 aircraft, a violent control surface buzz occurred and resulted in permanent damage to the control surface. Bendiksen (1993) noted other possible control surface buzz phenomena, which happened on Bell X-1 and X-1A aircrafts. The Bell X-1 aircraft encountered severe buffeting and control surface (aileron) buzz around Mach 0.86, however, as the aircraft accelerated over about Mach 0.95, the buzz phenomenon disappeared. A similar situation happened to X-1A aircraft, and the pilot saw the shocks rippling across the wing surface and felt the control

n

Corresponding author. E-mail addresses: [email protected] (S. He), [email protected] (Z. Yang), [email protected] (Y. Gu).

http://dx.doi.org/10.1016/j.jfluidstructs.2015.11.014 0889-9746/& 2016 Elsevier Ltd. All rights reserved.

432

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

surface (aileron) buzzing like mad. From above flight tests, we can see that the control surface buzz is dangerous, since it usually causes control surface vibrating suddenly and violently, and may result in a permanent damage on the structure. In the next two decades after the first observation of control surface buzz, wind tunnel tests were conducted to explore the physical mechanisms of this dangerous vibration. The study presented by Lambourne (1964) is one of the most important researches. Based on systematic wind tunnel tests, Lambourne suggested that control surface buzz could be classified into three types, namely Type A, Type B and Type C, according to Mach number and shock wave location on the control surface. Type A buzz typically occurs slightly above the critical Mach number when the shock stands forward of the control surface hinge, and this type buzz is thought to be related to shock-induced airflow separation. Type B buzz occurs at a higher Mach number, and it is associated with the shock oscillation on control surface. Type C buzz occurs at a even higher Mach number. The supersonic flow extends over the entire control surface and the shock attaches to the trailing edge. This classification for control surface buzz is still used in aircraft engineering nowadays. In 1990s,a series of transonic control surface buzz tests for National Aerospace Plane (NASP) wings were conducted in NASA Langley's Transonic Dynamic Tunnel. Parker et al. (1991) reported that control surface buzz is associated not only with Mach number but also with dynamic pressure, which disagrees with Lambourne's opinion that the occurrence of control surface buzz depends primarily on the Mach number rather than on the airspeed. Some of the experiment phenomena reported by Parker are of great interest. Some experiment models showed constant amplitude oscillations and divergent oscillations, and the other exhibited relatively low amplitude and heavily damped control surface oscillation just prior to divergent oscillations. Besides, the control surface buzz frequency increased by about 20–30% compared with a wind-off condition, which is called aerodynamic stiffening by Parker et al. (1991). Reasonable explanations for these interesting phenomena could contribute to understand the mechanisms of transonic control surface buzz. As a high-fidelity technique to simulating shock wave and flow separation in transonic flow, computational fluid dynamics (CFD) is usually employed to predict the control surface buzz. Steger and Bailey (1980) firstly used the numerical method to simulate control surface buzz based on Navier–Stokes equations, and their calculations indicated that control surface buzz did happen at a certain Mach number and showed good agreement with the results of wind tunnel tests. Pak and Baker (2001) applied CFL3D and CAPTSD codes to predict the initiation of control surface buzz at some Mach numbers according to those tested in wind tunnel. Constant amplitude buzz boundaries were captured during their simulation, and the buzz dynamic pressures and frequencies showed reasonably good agreement with those obtained by wind tunnel tests. A Hopf-bifurcation analysis method was employed to predict the transonic control surface buzz for a delta-wing model with a full-span control surface base on Navier–Stokes equations by Liu and Bai (2002), and the results achieved by the Hopfbifurcation analysis were consistent with those obtained by the time integration calculations and wind tunnel tests. A fully implicit multi-block aeroelastic solver, which coupled thin-layer Navier–Stokes equations with structural equations of motion, was implemented to simulate control surface buzz for a supersonic transport model by Yang and Obayashi (2003). Although the first six structural modes were taken into account in their simulations, the dominant mode was still the control surface oscillation mode when the control surface buzz occurred. With respect to Type B and Type C buzz, though boundary layer separation was also observed in the experiments, shock wave would be the dominant factor rather than separation (Bendiksen, 1993; Fusi et al., 2012). So solutions of Euler equations, which neglect the viscosity effect, are also sufficient to predict Type B and Type C control surface buzz. Bendiksen (1993) used Euler equations to simulate the aeroelastic behavior of a wing/control surface model, and the results suggested that control surface buzz could be caused by the interaction between shock motion and control surface vibration. In addition, wing flexibility in both bending and torsion, as well as camber bending, was modeled in his work. Based on Euler equations, Fusi et al. (2012) conducted a numerical study for Type B buzz aimed at drawing attention on grid details for accurate results. It can be seen from the above studies that Euler equations are also effective to simulate the control surface buzz of Type B and Type C. In essence, control surface buzz is a kind of single degree of freedom nonlinear flutter. According to Thomas et al. (2004), three possible types of limit cycle oscillation (LCO) behavior can be observed when nonlinearity is considered in the aeroelastic system, as illustrated in Fig. 1. The supercritical bifurcation is sometimes referred as soft flutter or benign LCO, in

Fig. 1. LCO behavior trends for nonlinear aeroelastic system from Thomas et al. (2004).

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

433

which no LCO occurs below the linear flutter speed. And the subcritical bifurcation is also called hard flutter or explosive LCO, where the LCO behavior trend consists of stable branch and unstable branch. The explosive LCO is dangerous because a small disturbance may be capable of triggering LCO even below the linear critical flutter speed. As we know, the unstable LCO cannot be obtained by wind tunnel test and time marching simulation, while frequency domain flutter solution can provide an entire subcritical bifurcation including both stable and unstable LCO branches according to Yang and Zhao (1988). Among the aforementioned numerical studies, the linear structural model was considered, and only the Mach number of control surface buzz onset was predicted by using CFD technique, but few researches on LCO behavior of transonic control surface buzz were conducted. Considering the free-play nonlinearity in the control surface hinge line is inevitable in a real aircraft, we also pay attention to the LCO behavior of transonic control surface buzz with both linear and nonlinear structural models in the current study. Aerodynamic describing function (He et al., 2014) and structural describing function (Yang and Zhao, 1988) are used to deal with aerodynamic and structural nonlinearities, respectively. And frequency domain flutter solutions are applied to obtain the stable and unstable LCOs. To the best of authors' knowledge, only Lambourne (1968) mentioned the nonlinear feature of transonic control surface buzz. The derivatives of the aerodynamic forces were measured in the wind tunnel by using forced-oscillation techniques in his work. He found that the aerodynamic damping was related to the oscillation amplitude and Mach number, which may lead to stable or unstable LCOs. Additionally, three typical examples of the variation of LCO amplitude with Mach number were shown, which can be regarded as the subcritical and supercritical bifurcations. However, the test model and relevant structure parameters are not available, and exact values of LCO amplitude and Mach number in bifurcation diagram are not illustrated in Lambourne (1968). Therefore, it is necessary to study the nonlinear feature of transonic control surface buzz in detail. In this paper, we demonstrate the LCO behavior of transonic control surface buzz. For physical simplicity, two simple but typical models, i.e. an all-movable control surface model and a wing/control surface model, are firstly adopted in the investigations. The allmovable control surface model is actually an isolated airfoil model with only pitching degree of freedom. It should be noted that the deflection angle of all-movable control surface is the pitching motion, but we also call it control surface deflection for convenience. A three-dimensional wing with a full-span control surface is also adopted to extend the demonstration. In the present study, Euler equations are used for Type B and Type C buzz analyses. Results for the control surface buzz system with linear structure show that the LCO behaviors of transonic buzz are subcritical or supercritical bifurcations at different Mach numbers. And the buzz frequency indicates that the aerodynamics has a stiffening effect on the aeroelastic system. For the nonlinear structural model with free-play in control surface deflection stiffness, the interesting double LCO phenomenon is observed. Comparing the results of linear and nonlinear structural models, it is found that the structural nonlinearity effect turns the LCOs at small amplitudes from unstable into stable. Moreover, the LCO behavior is dominated by the aerodynamic nonlinearity for large amplitude oscillation of control surface while by structural nonlinearity for small amplitude case.

2. Technical analysis As mentioned in Section 1, the transonic control surface buzz can be treated as a kind of single degree of freedom flutter, the solution of which can certainly be used to obtain the buzz speed and buzz frequency. 2.1. Single degree of freedom aeroelastic model The structural model used to simulate transonic buzz can only consider the deflection angle of control surface (Steger and Bailey, 1980; Fusi et al., 2012) or several structural modes including the mode of control surface deflection (Bendiksen, 1993; Yang and Obayashi, 2003). It is of importance that the results of Yang and Obayashi (2003) showed that the dominant mode was the control surface oscillation mode when the control surface buzz occurred, though the first six structural modes were taken into account. So in the current study, only the deflection angle of control surface β is considered, as sketched in Fig. 2. And the governing equation of motion can be written as I β β€ þ MðβÞ ¼ 2ρV 2 b cm ; 2

Fig. 2. Two-dimensional single degree of freedom flutter system in transonic flow.

ð1Þ

434

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449 2

where I β ¼ mr 2β b is the control surface moment of inertia about the hinge line, ρ is the air density, V is the air speed, b is the reference length which is usually the wing semi-chord length, cm is the aerodynamic moment coefficient about the hinge line, and m is the mass per unit span. Here β and cm are taken to be positive for downward deflection of the control surface. As usual, the hinge line is located at a distance of ab after midchord. In Eq. (1), MðβÞ represents the nonlinear structural restoring moment. When MðβÞ ¼ K β β, where K β ¼ I β ω2β is the control surface deflection stiffness about the hinge line, the problem is reduced to a linear structural model. But in the real aircraft wing or wind tunnel test model, free-play nonlinearity may exist in the control surface hinge line, which can be expressed as 8 K ðβ  δÞ; β Zδ; > < β δ oβ o δ; ð2Þ MðβÞ ¼ 0; > : K ðβ þ δÞ; β r δ; β where δ denotes the free-play region. 2 Introducing the mass ratio, μ ¼ m=πρb , Eq. (1) can be written as 1 € MðβÞ 2U 2 ¼ cm ; βþ Kβ ω2β πμr 2β

ð3Þ

where U ¼ V=bωβ is called the speed index. For this single degree of freedom aeroelastic system, β and cm are served as the generalized displacement and the generalized aerodynamic force, respectively. 2.2. Aerodynamic describing function Based on our previous work (He et al., 2014), aerodynamic describing function will be used to model the nonlinear aerodynamics for the present transonic single degree of freedom flutter analysis. As mentioned in Section 1, Euler equations are sufficient to simulate the Type B and Type C transonic buzz, and is also adopted to calculate the nonlinear aerodynamic force here. In the aerodynamic describing function method, Euler equations are regarded as an implicit nonlinear aerodynamic system between the generalized displacements (inputs) and the generalized aerodynamic forces (outputs). And the describing function is the ratio of the first order component of the complex output to the harmonic input with different amplitudes. The control surface deflection angle (input) is taken as a sinusoid function: βðtÞ ¼ β0 sin ðωtÞ;

ð4Þ

where β0 denotes the amplitude of control surface oscillation or LCO amplitude, and ω is the frequency. And the corresponding output of the CFD solver, aerodynamic moment coefficient, is periodic which can be expanded by Fourier series. In the present method, only the first order harmonic component of the aerodynamic moment coefficient is retained, and it can be written as cm ðtÞ ¼ ðcm Þc cos ðωtÞ þ ðcm Þs sin ðωtÞ ¼ ðcm Þcs sin ðωt þφÞ:

ð5Þ

The coefficient of ðcm Þc and ðcm Þs can be obtained by fitting the time history of aerodynamic moment to the above equation. Thus, the aerodynamic describing function can be expressed as Q¼

ðcm Þcs iφ ðcm Þs ðcm Þc e ¼ þi ¼ RðQ Þ þ iIðQ Þ: β0 β0 β0

ð6Þ

and the aerodynamic moment coefficient corresponding to the oscillation of control surface can be obtained by cm ¼ Q β:

ð7Þ

The aerodynamic describing function, Q, is related to Mach number, reduced frequency k ¼ ωb=V and the amplitude of control surface oscillation β0. It is obvious that the aerodynamic describing function is actually a kind of equivalent linearized aerodynamic influence coefficient corresponding to the amplitude of motion. 2.3. Structural describing function In our previous study (Yang and Zhao, 1988), harmonic balance analysis, which can also be called structural describing function method, was used to examine the LCO behavior of a two-dimensional wing model with nonlinear pitching stiffness. A fundamental harmonic solution for the control surface deflection angle is assumed, i.e. Eq. (4), then from Eq. (2) and the structural describing function approach, the linearized equivalent control surface deflection stiffness can be expressed as K eq ¼ f eq K β :

ð8Þ

In this expression, feq is called the structural describing function, which depends on the amplitude of control surface oscillation. And according to Gelb and Vander Velde (1968), the structural describing function for the free-play nonlinearity

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

is determined by 8 0; > > > 2 3 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <    2ffi 24 δ δ δ f eq ¼ 5; arcsin þ 1 1 > > > π β0 β0 β0 :

435

β0 o δ; ð9Þ

β0 Z δ:

Thus, the nonlinear structural restoring moment corresponding to the oscillation of control surface can be obtained by MðβÞ ¼ f eq K β β:

ð10Þ

2.4. Nonlinear flutter solution For nonlinear flutter problems, as Shen (1959) mentioned the aeroelastic system has such a strong filtering property that the fundamental harmonic dominates the flutter oscillations. Therefore, it is feasible to treat a simple harmonic oscillation in critical flutter condition: β ¼ β0 eiωt ;

ð11Þ

where β0 denotes the LCO amplitude, which is consistent with the definition in Eq. (4). Substituting Eqs. (7), (10) and (11) into the equation of motion Eq. (3) yields the flutter equation in frequency domain  2    ω 2U 2  β þf eq β0 β ¼ Q Ma; k; β0 β: ð12Þ  ωβ πμr 2β It is important to note that the structural describing function feq is a function of the LCO amplitude β0, and the aerodynamic describing function Q is a function of Mach number, reduced frequency k and β0. Thus, at a specific Mach number, the solution to Eq. (12) is not that of a standard linear flutter problem since β0 is unknown. An effective procedure to obtain the LCO of aeroelastic system with multiple lumped nonlinearities using describing function can be found in Manetti et al. (2009). However, considering that there is only one degree of freedom in the present model, a conventional procedure developed by Yang and Zhao (1988) and Shen (1959) is also applicable. For each prescribed LCO amplitude and at a specific Mach number, feq is a constant and Q is only related to reduced frequency k, so the nonlinear aeroelastic system is reduced to an equivalent linear aeroelastic system. For such a linear system, the flutter speed and flutter frequency can be obtained straightforwardly from the routine engineering flutter solution, for instance, V-g method or p-k method (Hodges and Pierce, 2011). The conventional V-g method is used to solve the flutter equation in the present study, and the artificial structural damping term D ¼ igβ is added to Eq. (12), yielding  2    ω 2U 2   β þigβ þ f eq β0 β ¼ Q Ma; k; β0 β: ð13Þ ωβ πμr 2β According to the procedure of V-g method, we obtain g¼

f eq ðβ0 ÞI½Q ðMa; k; β0 Þ 2

R½Q ðMa; k; β0 Þ þ πμr 2β k =2

;

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u f eq ðβ0 Þ ω=ωβ ¼ t ; 2 2R½Q ðMa; k; β0 Þ=πμr 2β k þ 1 U¼

ω=ωβ : k

ð14aÞ

ð14bÞ

ð14cÞ

Where ω=ωβ is called the frequency ratio. It is obvious that Eq. (14a) is an algebraic expression for a prescribed LCO amplitude at a specific Mach number. When the reduced frequency k is assigned, the corresponding g, ω=ωβ and U can be obtained. In the engineer applications of flutter analysis, g and ω=ωβ are usually plotted versus U to determine the critical flutter speed. According to Hodges and Pierce (2011), the plot of damping coefficient g versus U can indicate the margin of stability at conditions near the critical flutter speed. For a specific reduced frequency, the system is stable if damping coefficient is negative, and the system is unstable when damping coefficient is positive. When the damping coefficient approaches zero, the critical reduced frequency is achieved. And then the critical flutter frequency ratio and the critical flutter speed index can be obtained. Actually, the artificial structural damping g does not exist, and it is introduced to produce the desired harmonic motion. It is easy to find that when using conventional frequency domain flutter solution according to Hodges and Pierce (2011), such as V-g method, to solve the transonic flutter equation considering nonlinear aerodynamics, the only requirement is to replace the aerodynamic influence coefficient by the aerodynamic describing function. It should be pointed out that, by

436

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

setting f eq ¼ 1, the flutter solution (14a) is also suitable for transonic single degree of freedom flutter analysis with linear structure. By solving unsteady Euler equations and using describing function approach, the procedure of transonic aerodynamic nonlinear flutter analysis for a prescribed Mach number and a specific amplitude of control surface oscillation as defined in Eq. (4) can be described as follows: 1. 2. 3. 4. 5.

Obtain the structural describing function feq. Use a sinusoidal motion of the specific amplitude as the input of CFD solver for a prescribed reduced frequency. Solve unsteady Euler equations to calculate time history of the aerodynamic moment coefficient. Get the first order harmonic components by curve fitting of the computed aerodynamic moment coefficient. For several reduced frequencies, execute steps 2–4 to build the aerodynamic describing function, and then obtain the critical reduced frequency, critical flutter speed index and critical flutter frequency ratio.

Repeating the above procedure for different control surface oscillation amplitudes, the bifurcation diagram for transonic single degree of freedom flutter will be obtained. Conventionally, in the transonic control surface buzz analysis, the flutter speed index, flutter frequency ratio and the amplitude of the control surface is called as the buzz speed index, buzz frequency ratio and the buzz amplitude, respectively. 2.5. Stability of limit cycle oscillation The study of the stability of LCO is of great importance because the unstable LCO cannot be obtained by wind tunnel test and time marching simulation as mentioned in Section 1. A simple but effective way to obtain the stability of the LCO is to observe the curve of LCO amplitude versus flutter speed. According to Thomas et al. (2002) and Kholodar et al. (2004), the flutter speed decreases as the LCO amplitude increases, which suggests that this will lead to an unstable LCO, on the contrary the flutter speed increases as the LCO amplitude increases, which is indicative of a stable LCO. An alternative technique to determine the stability of the LCO using describing function is developed by Gelb and Vander Velde (1968) and Manetti et al. (2009). In frequency domain the flutter equation, Eq. (12), can be rewritten as "   #   2U 2   ω 2  þf eq β0  Q Ma; k; β0 β ¼ 0: ð15Þ ωβ πμr 2β From Eq. (14c), we obtain ω=ωβ ¼ kU:

ð16Þ

Taking Eq. (16) into Eq. (15) yields   2U 2   2 Q Ma; k; β0 ¼ 0: k U 2 þf eq β0  πμr 2β

ð17Þ

The above equation can be briefly indicated as Dðk; β0 Þ ¼ 0. According to Manetti et al. (2009), the LCO will be stable if it satisfies the following requirement: ∂R½D ∂I½D ∂R½D ∂I½D o0:  ∂k ∂β0 ∂β0 ∂k

ð18Þ

Note that the derivatives of Q, which is included in the derivatives of D, can be obtained numerically by using CFD solver.

3. LCO behaviors for linear structural model To verify the accuracy of the CFD solver used in the present study, unsteady aerodynamic forces acting on an oscillating NACA 0012 airfoil are computed and compared to the experimental results of AGARD CT Case 5 (Landon, 1982). The computed Mach number is 0.755, and this simulation considers the prescribed oscillating motion of the airfoil given in the following equation: βðtÞ ¼ βm þ β0 sin ðωtÞ;

ð19Þ

where βm ¼ 0:0161, β0 ¼ 2:511, ω ¼ 36:278 rad=s. The rotation axis is located at 1/4 chord point of the airfoil. The computed lift and moment coefficient versus pitching angle and the experimental data are depicted in Fig. 3, in which good agreement is shown. In this section, the linear structural model is considered. Firstly, we select an all-movable control surface model to verify the accuracy of the present transonic single degree of freedom flutter solution. Then, we investigate the LCO behaviors of Type B and Type C control surface buzz for the all-movable control surface model and the wing/control surface model. As described in Section 1, the all-movable control surface model is de facto an isolated airfoil model with only pitching degree of freedom.

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449 0.8

0.05

0.6

0.04 0.03

0.4

cl

0.01 0

cm

0.02

0.2

0

-0.2

-0.01

-0.4 -0.6

437

-0.02 -3

-2

-1

0

1

2

3

-0.03

β (deg) Fig. 3. Comparison of computed lift/moment coefficient with experimental data.

Fig. 4. Computational grids for NACA 64A010: (a) overall and (b) close-up.

3.1. Test case The all-movable control surface model of the test case is taken from the work of Thomas et al. (2002), in which the harmonic balance (HB) method is used to calculate the transonic unsteady aerodynamic force. The airfoil of this model is NACA 64A010, and its parameters are μr 2β ¼ 50; a ¼ 0:6. The computed Mach number is 0.8. Three computational grids with different spatial resolutions (Mesh 1, Mesh 2 and Mesh 3) are generated to assess the mesh convergence, where the C-type structured mesh topology is adopted. All these three computational meshes are composed of 51 mesh points radially and the outer boundary for computational domain extends to a distance of 50 chord length from the airfoil, as shown in Fig. 4. However, the number of the mesh points surrounding the airfoil surface is different, for instance, 111 mesh points in Mesh 1, 221 in Mesh 2 and 331 in Mesh 3. So there are totally 211, 321 and 431 mesh points circumferentially in Mesh 1, Mesh 2 and Mesh 3, respectively. CFD simulations of the all-movable control surface model with a sinusoid oscillating motion for β0 ¼ 11 and k ¼ 0:1 at Mach 0.8 are conducted. The time histories of the unsteady aerodynamic coefficient for different meshes are shown in Fig. 5. It can be seen that the amplitude of cm using Mesh 1 is a little larger than those using Mesh 2 and Mesh 3, but results obtained by using Mesh 2 and Mesh 3 are nearly the same. It reveals that both Mesh 2 and Mesh 3 are feasible to calculate the unsteady aerodynamic forces. Considering the computational costs, Mesh 2 is adopted in this section. At a specific amplitude of control surface oscillation, after building the aerodynamic describing function for a series of reduced frequencies, V-g method is used to calculate the critical flutter speed index and the flutter frequency ratio. Depicting these results versus the amplitude of control surface deflection angle, the bifurcation diagram is obtained, as illustrated in Fig. 6(a) and (b). As compared with the bifurcation diagram obtained by HB method from Thomas et al. (2002), the LCO behavior trends are in reasonably good agreement, which demonstrates the feasibility of using the present method for transonic LCO solution. Thomas et al. (2002) noted that the curve of LCO amplitude versus speed index bending to the left indicates the unstable LCO branch, the stability of which agrees with that determined by Eq. (18). As mentioned in Section 1, the stable LCO branch for subcritical bifurcation should exist when the amplitude of control surface is further increased. To seek for the stable LCO branch, the flutter solutions for the cases with larger amplitude up to 18.2° are calculated, as shown in Fig. 6(c) and (d). As

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

0.03

Mesh 1 Mesh 2 Mesh 3

0.02

0.01

cm

438

0

-0.01

-0.02

-0.03

0

0.1

0.2

0.3

0.4

0.5

0.6

Time (s) Fig. 5. Comparison of aerodynamic moment coefficient of NACA 64A010 for different meshes (Mach 0.8, β0 ¼ 11, k ¼ 0.1 and a ¼  0:6).

Fig. 6. Comparison of the bifurcation diagram for NACA 64A010 (Mach 0.8, μr 2β ¼ 50): (a) and (c) speed index and (b) and (d) frequency ratio.

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

439

we expect, when the amplitude of control surface oscillation is increased greater than 16.5°, the curve of LCO amplitude versus speed index turns to bending right, implying the stable branch. So the LCO behavior for the test case is a type of subcritical bifurcation consisting of both stable and unstable branches. 3.2. All-movable control surface model The airfoil of this all-movable control surface model is also NACA 64A010. The computational grids of Mesh 2 in Section 3.1 are still used here. And the other parameters are μ ¼ 75; r 2β ¼ 0:75; a ¼  0:6, which are taken from Thomas et al. (2002) for transonic LCO analysis of a two-dimensional wing. Fig. 7 shows the steady pressure coefficient distribution for NACA 64A010 airfoil at control surface deflection angle of 0° but at Mach 0.875, 0.91 and 0.95. From the classification by Lambourne (1968), there might be only Type B and Type C buzz happening on the so-called all-movable control model, since the shock either lies on the airfoil or attaches to the trailing edge when it forms. This classification does not seem to be very meaningful for the all-movable control surface model, but it should be emphasized that shock movement is the dominant factor rather than flow separation at these three Mach numbers, so the adoption of Euler equations is adequate in the present study. Fig. 8 shows the bifurcation diagram at Mach 0.875 and results from time marching simulation based on Euler equations. There are good agreements between results of the frequency domain solution and the time marching approach, demonstrating the validity of aerodynamic describing function for buzz analysis. As mentioned in Section 1, the time marching approach can only simulate the stable LCOs but cannot obtain the unstable LCOs, which can be illustrated in Fig. 8. -1

-0.5

cp

0

0.5

Mach 0.85 Mach 0.875 Mach 0.9 Mach 0.91 Mach 0.95

1

1.5

0

0.2

0.4

0.6

0.8

1

x Fig. 7. Steady pressure coefficient distribution at different Mach numbers (β ¼ 01).

Fig. 8. Bifurcation diagram for the all-movable control surface model at Mach 0.875: (a) speed index and (b) frequency ratio. Solid line and dashed line stand for stable and unstable branches, respectively. Cycle (◯) stands for results from time marching method based on CFD.

440

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

Fig. 9 shows the bifurcation diagram at Mach 0.91. It can be seen that the buzz speed increases with the increasing of LCO amplitude, which indicates a supercritical bifurcation. The buzz frequency ratio ω=ωβ can reflect the effect of aerodynamics on the aeroelastic system. If the buzz frequency ratio is greater than 1, the aerodynamics has a stiffening effect, i.e. aerodynamic stiffening effect, on the aeroelastic system (Parker et al., 1991). The buzz frequency ratio at this Mach number is greater than 1.04, implying the aerodynamic stiffening effect. Fig. 10 shows the bifurcation diagram at Mach 0.95. As shown in Fig. 10(b), the critical buzz frequency for stable branch is in the range from 1.2 to 1.3, which is consistent with wind tunnel experiment observation that the control surface frequency increased by about 20–30% compared with a wind-off condition according to Parker et al. (1991). From Figs. 8–10, it is observed that the LCO behavior for the all-movable control surface model can be either subcritical bifurcation (at Mach 0.875 and 0.95) or supercritical bifurcation (at Mach 0.91). For the type of supercritical bifurcation, stable LCOs can oscillate with a large range of amplitude from 1° to 18°. But for the type of subcritical bifurcation, the amplitude of stable LCO is in the narrow range from 17° to 23°. 3.3. Wing/control surface model In this section, a wing/control surface model is analyzed, whose airfoil is NACA 0012. Other parameters are μ ¼ 25:24; r 2β ¼ 0:01292; a ¼ 0:5. The structural parameters are taken from Dowell et al. (2004) for transonic LCO analysis considering both structural and aerodynamic nonlinearities of a two-dimensional wing/control surface system.

Fig. 9. Bifurcation diagram for the all-movable control surface model at Mach 0.91: (a) speed index and (b) frequency ratio. Solid line stands for stable branch.

Fig. 10. Bifurcation diagram for the all-movable control surface model at Mach 0.95: (a) speed index and (b) frequency ratio. Solid line and dashed line stand for stable and unstable branches, respectively.

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

441

Similar to Section 3.1, the C-type structured mesh topology is used and four computational grids (Mesh 1–Mesh 4) with different spatial resolutions are generated to evaluate the mesh convergence. The number of the mesh points on the control surface is different, for instance, 28 points in Mesh 1, 56 points in Mesh 2, 84 points in Mesh 3 and 112 points in Mesh 4. So there are totally 293, 321, 349 and 377 grid points circumferentially in these four meshes, respectively. The computational grids are composed of 51 mesh points radially and the outer boundary for computational domain extends to a distance of 50 chord length from the airfoil for all the four meshes, as shown in Fig. 11. CFD simulations of the wing/control surface model with a sinusoid oscillating motion for β0 ¼ 11 and k ¼ 0:1 at Mach 0.86 are implemented. The time histories of the unsteady aerodynamic coefficient for different meshes are compared in Fig. 12. It can be seen that as the number of grid points on the control surface increases from Mesh 1 to Mesh 4, the difference of the time history of cm becomes smaller. Considering the computational costs, Mesh 3 is adopted in this section. Fig. 13 shows the steady Mach number contour for the wing/control surface model at Mach 0.86, 0.93 and 1.2 with control surface deflection angle of 2°. It can be seen that the shock wave forms on the control surface at Mach 0.86, and the supersonic region extends over the entire control surface and the shock wave occurs at the control surface trailing edge at Mach 0.93 and 1.2. Figs. 14–16 show the bifurcation diagram at Mach 0.86, 0.93 and 1.2, respectively. Similar to the all-movable control surface model, the subcritical bifurcation and supercritical bifurcation are observed at different Mach numbers. In addition, Figs. 15 and 16 show the LCOs from the time marching approach based on Euler equation, where good agreements between the present method and the time marching method are illustrated. It should be noted that no LCO is detected if the LCO amplitude is less than 10° at Mach 1.2 as shown in Fig. 16. To sum up, from all the results of the linear structural model, it is found that the bifurcation of the transonic control surface buzz can be subcritical or supercritical at different Mach numbers for the all-movable control surface model and the

Fig. 11. Computational grids for wing/control surface model: (a) overall and (b) close-up.

Mesh 1 Mesh 2 Mesh 3 Mesh 4

0.006

0.004

cm

0.002

0

-0.002

-0.004

-0.006

0

0.1

0.2

0.3

0.4

0.5

0.6

Time (s) Fig. 12. Comparison of aerodynamic moment coefficient of wing/control surface model for different meshes (Mach 0.86, β0 ¼ 11, k ¼ 0:1 and a ¼ 0:5).

442

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

Fig. 13. Steady Mach number contour for wing/control surface model at β ¼ 2○ : (a) Mach 0.85, (b) Mach 0.93, and (c) Mach 1.2.

Fig. 14. Bifurcation diagram for the wing/control surface model at Mach 0.86: (a) speed index and (b) frequency ratio. Solid line and dashed line stand for stable and unstable branches, respectively.

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

443

Fig. 15. Bifurcation diagram for the wing/control surface model at Mach 0.93: (a) speed index and (b) frequency ratio. Solid line stands for stable branch. Cycle (◯) stands for results from time marching method based on CFD.

Fig. 16. Bifurcation diagram for the wing/control surface model at Mach 1.2: (a) speed index and (b) frequency ratio. Solid line and dashed line stand for stable and unstable branches, respectively. Cycle (◯) stands for results from time marching method based on CFD.

wing/control surface model. Let us recall the experiment observation that the model exhibited heavily damped control surface oscillation just before divergent oscillations reported by Parker et al. (1991). And with the knowledge of subcritical bifurcation for transonic control surface buzz in the present study, the experimental observation might be easy to understand. Moreover, the experiment observation that the model exhibited constant amplitude oscillations before divergent is also understandable, if the bifurcation type for the experiment model is supercritical.

4. LCO behaviors for nonlinear structural model 4.1. Bifurcation diagrams A free-play nonlinearity with δ ¼ 0:51 is assumed in the stiffness of the control surface deflection herein. Fig. 17 shows the bifurcation diagrams considering both aerodynamic and structural nonlinearities for the all-movable control surface model and the wing/control surface model at different Mach numbers. It can be seen that for the all-movable control surface model at Mach 0.91 and the wing/control surface model at Mach 0.93, types of bifurcation are still supercritical. However, for other two Mach numbers of each model, the segment of unstable branch with small amplitudes in linear structural model turns to be stable when the structural nonlinearity is taken into account. Thus, the stable LCOs with not only large amplitudes but

444

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

Fig. 17. Bifurcation diagram considering both aerodynamic and structural nonlinearities: (a), (c) and (e) all-movable control surface model at Mach 0.875, Mach 0.91 and Mach 0.95, respectively, and (b), (d) and (f) wing/control surface model at Mach 0.86, Mach 0.93 and Mach 1.2, respectively. Solid line and dashed line stand for stable and unstable branches, respectively.

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

445

also small amplitudes are observed. In the following section, we will pay more attention on how the cases of subcritical bifurcation will behave when the free-play nonlinearity is considered. As shown in Fig. 17(e), when speed index is less than U1, only LCOs with small amplitude will be observed. If the speed index is increased between U1 and U2, there are three LCO amplitudes, A1, A2 and A3, but only LCOs with small amplitude A1 and large amplitude A3 are stable, which is termed as double LCOs phenomenon. This interesting phenomenon was also observed in our previous study of Yang and Zhao (1988). From the authors' experience, with adequate initial condition in this speed index range, LCOs with two different amplitudes will be observed in the wind tunnel tests and time domain simulations. When the speed index is increased further, LCOs at large amplitudes are shown, the amplitudes of which are very close to those of linear structure. For a real aircraft, the double LCOs phenomenon could be very dangerous, because the LCO with small oscillation amplitudes may turn to vibrate with very large amplitude with perturbations, which may result in a permanent structural damage. We can imagine that if the free-play nonlinearity exists in the control surface hinge line of a real aircraft or wind tunnel test model, as the air dynamic pressure increases, LCOs with large amplitudes could appear abruptly after the emergence of the LCOs with small amplitudes. This phenomenon reminds us of the experiment observation that some model behaves relatively low amplitude control surface oscillation just before divergent oscillations reported by Parker et al. (1991). Our research for transonic control surface buzz considering both aerodynamic and structural nonlinearities might be a possible physical mechanism explaining the experiment phenomenon. 4.2. Discussion For the cases of all-moveable model (Mach 0.875 and 0.95) and the wing/control surface model (Mach 0.86 and 1.2), the bifurcation diagrams between the case considering only aerodynamic nonlinearity and the case considering both aerodynamic and structural nonlinearities are compared, and it can be found that the structural nonlinearity changes the LCOs with small amplitude from unstable into stable. To address this interesting phenomenon, the relationship between buzz speed index of linear structure Ul and buzz speed index for nonlinear structure Unon is explored from the V-g solution. From the flutter solution mentioned in Section 2.4, it can be seen that we use the condition g ¼ 0 to determine the critical reduced frequency for a prescribed buzz amplitude. And from Eq. (14a), we know that the condition g ¼ 0 only depends on IðQ Þ ¼ 0, which is determined by the aerodynamic model. So the structural parameters do not have an influence on the critical reduced frequency. That is to say, the critical reduced frequencies at a specific buzz amplitude are the same for both linear and nonlinear structures. Substituting this critical reduced frequency kcr and corresponding aerodynamic moment coefficient Qcr into Eqs. (14b) and (14c) yields vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u qffiffiffiffiffiffi f eq ðω=ωβ Þcr u U non ¼ =kcr ¼ f eq U l : ¼t ð20Þ 2 2 kcr 2RðQ cr Þ=πμr k þ1 β cr

Thus, the buzz speed index for nonlinear structure Unon can be obtained by scaling the buzz speed index for linear structure qffiffiffiffiffiffi Ul with the factor f eq at the same buzz amplitude. With respect to free-play nonlinearity considered in the present study, Fig. 18 shows the relationship between scale qffiffiffiffiffiffi factor f eq and amplitude of control surface deflection angle. It can be seen that the scale factor values are less than 1 for all the oscillation amplitudes of control surface, which reveals that the buzz speed index considering free-play nonlinearity is smaller than that for linear structural model at the same amplitude. As the amplitude increases, since the scale factor value

Fig. 18. Scale factor

qffiffiffiffiffiffi f eq versus amplitude of control surface deflection angle (δ ¼ 0.5 1).

446

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

Fig. 19. Comparison of LCO behaviors between with and without free-play nonlinearity (all-movable control surface model at Mach 0.95).

Fig. 20. Description of NASP wing model with a full-span control surface: (a) baseline model and (b) model A.

increases, the buzz speed will also increase if only the free-play nonlinearity is considered, resulting in a stable LCO branch. When the amplitude is increased larger than 4°, however, the scale factor value approaches nearly 0.95, which indicates that the effect of structural nonlinearity becomes very weak at this stage. As shown in Fig. 19, when the amplitude of control surface motion is less than Acr, the buzz speed for linear structural model decreases with the increasing of amplitude. But if the free-play nonlinearity is considered in the structural model, the change of the buzz speed with the increasing of control surface oscillation amplitude becomes complex and relies on the amplitude. As mentioned above, when the amplitude is larger than 4° in the present study, the structural nonlinearity is weak, and the LCO behavior is dominated by the aerodynamic nonlinearity. Consequently, the bifurcation diagram behaves like control surface buzz of the linear structure, and the stable LCO branch with large amplitudes and unstable LCOs branch are observed. On the contrary, if the amplitude is less than 4°, the structural nonlinearity dominates the LCO behavior, resulting in the stable LCOs with small amplitudes.

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

447

Table 1 Structural parameter values for NASP baseline model. Dimensional quantities

Non-dimensional quantities

For full span length

For unit span length

m ¼ 0:267788 kg I β ¼0.000881959 kg m2

m ¼ 0:857141 kg m  1 Iβ ¼0.00282299 kg m

μ ¼ 1:53432 r2β ¼0.0226886 a ¼ 0:726667

b ¼ 0:381 m ωβ ¼ 29:0 Hz

Table 2 Structural parameter values for NASP model A. Dimensional quantities

Non-dimensional quantities

For full span length

For unit span length

m ¼ 0:190629 kg Iβ ¼ 0.000327430 kg m2

m ¼ 0:641458 kg m  1 Iβ ¼ 0.00110179 kg m

b ¼ 0:508 m ωβ ¼ 31:6 Hz

μ¼ 1.58242 r2β ¼ 0.00665584 a ¼ 0:8075

5. LCO behaviors for NASP model The LCO behaviors of the baseline model and model A of NASP wings with a full-span aileron according to Parker et al. (1991) are examined in this section. The planforms of these two models are shown in Fig. 20, and the airfoil for both models is a 3% thick biconvex airfoil. The semi-chord length at the wing root is taken as the reference length, and the aileron hinge line is vertical to the wing root chord. Other relevant structural parameters for the aileron are listed in Tables 1 and 2. Since the mass mn and inertial moment I β mentioned in Parker et al. (1991) are measured for full span length, we transform them to m and Iβ for unit span by dividing the length of span, and then to non-dimensional parameters which are suitable for buzz analysis in the present study. Fig. 21 shows the bifurcation diagram for baseline model at Mach 1 and the experimental results from Parker et al. (1991). Since no values of LCO amplitude are available in Parker et al. (1991), two lines parallel to the axis of the LCO amplitude are used to represent the buzz dynamic pressure for comparison. The experimental observation of transonic buzz is described as that the baseline model exhibited constant amplitude oscillations prior to divergent oscillations by Parker et al. (1991). Our numerical results show that the type of bifurcation for baseline model at Mach 1 is supercritical. With the knowledge of LCO behavior for the baseline model, the constant amplitude buzz and divergent amplitude buzz may be easy to understand. Parker et al. (1991) mentioned that the nature of the oscillation on model A did vary from those observed on the baseline model. Specially, model A exhibited relatively low amplitude and heavily damped aileron oscillations just prior to divergent oscillations. Based on the previous discussions, it can be inferred that LCO behavior for model A might be subcritical bifurcation. Fig. 22 shows the bifurcation diagram for model A at Mach 0.97 from the present method and the dynamic pressure range in which divergent amplitude buzz was observed by Parker et al. (1991). As we expect, the type of bifurcation for model A is subcritical bifurcation. We can imagine that if the structural nonlinearity, such as free-play, exists in the control surface deflection stiffness, then LCOs with small amplitudes may be observed in the experiment. Comparing wind tunnel experimental observations and LCO behaviors obtained by the present method, it can be seen that our study for transonic buzz may be helpful to explain the experimental observations in wind tunnel tests.

6. Conclusions Euler equations are employed to calculate the unsteady aerodynamic forces for Type B and Type C transonic control surface buzz analyses, and an all-movable control surface model, a wing/control surface model and a three-dimensional wing with a full-span control surface are adopted in the present study. Aerodynamic describing function and structural describing function are used to deal with aerodynamic and structural nonlinearities, respectively. Then the buzz speed index and buzz frequency ratio are obtained by V-g method. For linear structural model, the LCO behavior exhibits subcritical bifurcation or supercritical at different Mach numbers. The subcritical bifurcation consists of both stable and unstable LCO branches, which suggests that the transient aeroelastic response of control surface for a small disturbance will be damped and then oscillate with large amplitudes as the air

448

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

Fig. 21. Bifurcation diagram for NASP baseline model at Mach 1: (a) dynamic pressure and (b) frequency. Solid line stands for stable branch.

Fig. 22. Bifurcation diagram for NASP model A at Mach 0.97: (a) dynamic pressure and (b) frequency. Solid line and dashed line stand for stable and unstable branches, respectively.

dynamic pressure increases. The supercritical bifurcation implies that with the increasing of air dynamic pressure the transient aeroelastic response of control surface will be constant amplitude oscillations and then gradually oscillate with large amplitudes. For nonlinear structural model with free-play nonlinearity in the control surface deflection stiffness, the double LCO phenomenon is observed in certain range of speed index. The bifurcation diagrams imply that as the air dynamic pressure increases, the control surface oscillation with large amplitudes may take place suddenly after the emergence of the LCOs with small amplitudes. Besides, when free-play nonlinearity is considered, the stability of LCOs at small amplitudes changes into stable from unstable for linear structure case. For large amplitudes of control surface oscillations, the LCO behavior is dominated by the aerodynamic nonlinearity, while the LCO behavior with small amplitudes is dominated by the structural nonlinearity. Our study for transonic buzz with both linear and nonlinear structural model may be helpful to explain the experimental observations in wind tunnel tests and to understand the physical mechanism of transonic control surface buzz.

Acknowledgements This work is supported by the National Natural Science Foundation of China (Grant nos. 11472216 and 11102162) and the 111 Project of China (Grant no. B07050).

S. He et al. / Journal of Fluids and Structures 61 (2016) 431–449

449

References Bendiksen, O.O., April 1993. Nonclassical aileron buzz in transonic flow. In: 34th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, No. AIAA-93-1479. American Institute of Aeronautics and Astronautics La Jolla, CA, USA. Dowell, E.H., Thomas, J.P., Hall, K.C., 2004. Transonic limit cycle oscillation analysis using reduced order aerodynamic models. Journal of Fluid and Structures 19, 17–27. Fusi, F., Romanelli, G., Guardone, A., Quaranta, G., November 2012. Assessment of geometry reconstruction techniques for the simulation of non-classical aileron buzz. In: ASME 2012 International Mechanical Engineering Congress and Exposition, No. IMECE 2012-87464. American Society of Mechanical Engineers, Houston, Texas, USA. Gelb, A., Vander Velde, W.E., 1968. Multiple-Input Describing Functions and Nonlinear System Design. McGraw-Hill Book Company, New York. He, S., Yang, Z., Gu, Y., 2014. Transonic limit cycle analysis using aerodynamic describing function and superposition principle. AIAA Journal 52, 1393–1403. Hodges, D.H., Pierce, G.A., 2011. Introduction to Structural Dynamics and Aeroelasticity. Cambridge University Press, New York. Kholodar, D.B., Dowell, E.H., Thomas, J.P., Hall, K.C., 2004. Limit-cycle oscillation of a typical airfoil in transonic flow. Journal of Aircraft 44 (5), 1067–1072. Lambourne, N.C., 1964. Control-surface buzz. In: Aeronautical Research Council Reports and Memoranda. No. 3364. Ministry of Aviation, Her Majestys' Stationery Office, London, UK. Lambourne, N.C., February 1968. Flutter in one degree of freedom. In: AGARD Manual on Aeroelasticity, vol. 5. Advisory Group for Aerospace Research & Development, France. Landon, R.H., August 1982. NACA0012 Oscillatory and Transient Pitching, No. AGARD-R-702. Liu, Q., Bai, J., September 2002. Transonic aileron buzz boundary prediction using Hopf-bifurcation analysis method. In: 23rd Congress of International Council of the Aeronautical Sciences, No. ICAS 2002-2.7.4. Toronto, Canada. Manetti, M., Quaranta, G., Mantegazza, P., 2009. Numerical evaluation of limit cycles of aeroelastic systems. Journal of Aircraft 46 (5), 1759–1769. Pak, C.-g., Baker, M.L., April 2001. Control surface buzz analysis of a generic NASP wing. In: 42nd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, No. AIAA-2001-1581. American Institute of Aeronautics and Astronautics, Seattle, WA, USA. Parker, E.C., Spain, C.V., Soistmann, D.L., 1991. Aileron buzz investigated on several generic nasp wing configurations. In: 32nd Structures, Structural Dynamics, and Materials Conference, No. AIAA-91-0936-CP. American Institute of Aeronautics and Astronautics, Baltimore, MD, USA. Shen, S.F., 1959. An approximate analysis of nonlinear flutter problems. Journal of the Aerospace Sciences 28, 25–32. Steger, J.L., Bailey, H.E., 1980. Calculation of transonic aileron buzz. AIAA Journal 18, 249–255. Thomas, J.P., Dowell, E.H., Hall, K.C., 2002. Nonlinear inviscid aerodynamic effects on transonic divergence, flutter and limit cycle oscillations. AIAA Journal 40, 638–646. Thomas, J.P., Dowell, E.H., Hall, K.C., 2004. Modeling viscous transonic limit-cycle oscillation behavior using a harmonic balance approach. Journal of Aircraft 41, 1266–1274. Yang, G., Obayashi, S., 2003. Aileron buzz simulation using an implicit multiblock aeroelastic solver. Journal of Aircraft 40, 580–589. Yang, Z.C., Zhao, L.C., 1988. Analysis of limit cycle flutter of an airfoil in incompressible flow. Journal of Sound and Vibration 123, 1–13.