Transonic limit cycle oscillation behavior of an aeroelastic airfoil with free-play

Journal of Fluids and Structures 66 (2016) 1–18

Contents lists available at ScienceDirect

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

Transonic limit cycle oscillation behavior of an aeroelastic airfoil with free-play Zhichun Yang n, Shun He, 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 28 October 2015 Received in revised form 4 July 2016 Accepted 19 July 2016

The limit cycle oscillation (LCO) behaviors of an aeroelastic airfoil with free-play for different Mach numbers are studied. Euler equations are adopted to obtain the unsteady aerodynamic forces. Aerodynamic and structural describing functions are employed to deal with aerodynamic and structural nonlinearities, respectively. Then the flutter speed and flutter frequency are obtained by V-g method. The LCO solutions for the aeroelastic airfoil obtained by using dynamically linear aerodynamics agree well with those obtained directly by using nonlinear aerodynamics. Subsequently, the dynamically linear aerodynamics is assumed, and results show that the LCOs behave variously in different Mach number ranges. A subcritical bifurcation, consisting of both stable and unstable branches, is firstly observed in subsonic and high subsonic regime. Then in a narrow Mach number range, the unstable LCOs with small amplitudes turn to be stable ones dominated by the single degree of freedom flutter. Meanwhile, these LCOs can persist down to very low flutter speeds. When the Mach number is increased further, the stable branch turns back to be unstable. To address the reason of the stability variation for different Mach numbers at small amplitude LCOs, we find that the Mach number freeze phenomenon provides a physics-based explanation and the phase reversal of the aerodynamic forces will trigger the single degree of freedom flutter in the narrow Mach number range between the low and high Mach numbers of the chimney region. The high Mach number can be predicted by the freeze Mach number, and the low one can be estimated by the Mach number at which the aerodynamic center of the airfoil lies near its elastic axis. Influence of angle of attack and viscous effects on the LCO behavior is also discussed. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Transonic flow Aerodynamic nonlinearity Free-play Limit cycle oscillation Transonic flutter

1. Introduction Limit cycle oscillation in aeroelastic system seems to be more prevalent in transonic air flow than in subsonic air flow according to Thomas et al. (2002). Certainly, LCO may be caused by structural nonlinearity no matter whether the air flow is transonic or not (Yang and Zhao, 1988), and the aerodynamic nonlinearity will lead to LCOs in transonic air flow even in the absence of any structural nonlinearity (Thomas et al., 2002). However, the research on a high aspect ratio swept wing conducted by Bendiksen (2008) showed that if a linear structural model is used in transonic flutter analysis, the critical dynamic pressure at LCO onset is over-predicted by a factor of about 3 compared with that observed in wind-tunnel tests. It revealed the importance of using a nonlinear structural model in transonic LCO analysis. In this paper, the LCO behaviors of n

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

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

2

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

an aeroelastic airfoil with free-play nonlinearity in transonic air flow are discussed. Usually, the aeroelastic airfoil with concentrated nonlinearities and the high aspect ratio wing with geometrical nonlinearity are adopted as the theoretical models studied in transonic air flow. Kousen and Bendiksen (1994) investigated the transonic LCO of an aeroelastic airfoil with a free-play nonlinearity in the pitching degree of freedom (DOF) and the bifurcation diagram was obtained by integrating the motion of the airfoil, i.e. the time marching approach. The transonic aeroelastic behaviors of an airfoil with free-play nonlinearities in both of the pitching and plunging DOFs were analyzed by Kim and Lee (2000), and LCOs and chaotic motions were observed in specific ranges of Mach numbers and air speeds. A three DOFs aeroelastic wing model with free-play in its control surface deflection was used to study the transonic LCO behavior by Dowell et al. (2004b), and the so-called chimney phenomenon was shown in their results. “Chimney” is a term used to describe the rapid change of the flutter dynamic pressure with the increasing of Mach number and it is often associated with a sudden change in flutter mode. The LCOs for an aeroelastic airfoil with structural nonlinearity in subsonic and transonic air flow were obtained by Munteanu et al. (2005) using the reduced order model (ROM) method, and their results showed that the ROM approach can estimate the LCOs accurately and efficiently. An airfoil with piecewise nonlinearity was investigated by Jones et al. (2007) in both low speed incompressible air flow and transonic air flow, and the LCOs were rapidly identified with good accuracy in their study. A delta wing model with geometrical nonlinearity in transonic air flow was studied by Chen et al. (2014), and a new nonlinear aeroelastic model was developed in their work. All of the research activities mentioned above focused on predicting the response of the nonlinear aeroelastic systems precisely. However, the LCO behaviors for different Mach numbers in transonic air flow were not demonstrated systematically. In this study, we use an aeroelastic airfoil with free-play in its pitching DOF to explore how the LCO behaves at different Mach numbers. According to Thomas et al. (2004), three possible types of LCO behavior can be observed in nonlinear aeroelastic wing systems, as illustrated in Fig. 1. The supercritical bifurcation is sometimes referred as the soft flutter or benign LCO, in which no LCO occurs below the linear critical flutter speed. And the subcritical bifurcation is also called the hard flutter or explosive LCO, where the LCO has a stable branch and an unstable branch. In Fig. 1, the curve of LCO amplitude versus airspeed for the so-called linear LCO behavior is perpendicular to the airspeed axis. It should be noted that the time marching approach can only capture the stable LCOs, but a frequency domain flutter solution can obtain both stable and unstable LCOs forming the entire LCO behaviors according to Yang and Zhao (1988). The time marching approach with computational fluid dynamics (CFD) technique (Kousen and Bendiksen, 1994; Bendiksen, 2008; Kim and Lee, 2000; Chen et al., 2014) or the ROM method based on CFD technique (Jones et al., 2007; Munteanu et al., 2005) is conventionally used to simulate the transonic LCOs. But as we know the unstable LCOs can not be obtained by using the time marching approach. Therefore, in the present study the aerodynamic describing function (He et al., 2014) and structural describing function (Yang and Zhao, 1988) are adopted to deal with the aerodynamic and structural nonlinearities, respectively. And the frequency domain flutter solution is subsequently applied to obtain the stable and unstable LCOs. When the shock wave is generated in transonic air flow, the steady air flow parameters vary with spatial position in the flow field around the wing. So the behaviors of the unsteady aerodynamic forces induced by the airfoil motion will be nonlinear in the transonic region. Nevertheless, it can be assumed that the parameters of flow field and the motion of shock wave vary in a linear fashion with the wing motion when considering a small perturbation about the transonic steady air flow. And this is usually called dynamically linear but statically nonlinear aerodynamics according to Dowell et al. (2004a). In the current study, Euler equations are used to calculate the unsteady aerodynamic forces for transonic LCO analysis. The dynamically linear aerodynamics is adopted to model the transonic aerodynamics for LCO analysis at different Mach numbers. Our results show that in subsonic and high subsonic air flows, the subcritical bifurcation is observed, consisting of

Fig. 1. LCO behaviors for nonlinear aeroelastic wing from Thomas et al. (2004).

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

3

Fig. 2. An aeroelastic airfoil in transonic air flow.

both stable and unstable branches. As the Mach number is increased, the unstable branch for small LCO amplitudes turns to be stable in a narrow Mach number range. The LCO can be observed at a very low flutter speed at this Mach number range. With the further increasing of Mach number, the stable branch turns back to be unstable. Our investigation here will focus on why the unstable LCOs with small amplitudes turn to be stable ones in this specific Mach number range and how to estimate the low and high Mach numbers of this range.

2. Technical analysis 2.1. Governing equations for an aeroelastic airfoil Fig. 2 shows a sketch of a typical aeroelastic airfoil with plunging and pitching DOFs, and the governing equation of motion can be written as

⎧ ¨ 2 ⎪ mh + S α α ¨ + Kh h = − ρV bcl ⎨ ⎪ 2 ¨ ⎩ Sα h + Iα α¨ + M (α ) = 2ρV b2cm

(1)

where h and α are the airfoil plunging and pitching, cl and cm are the lift coefficient and aerodynamic moment coefficient about the elastic axis, m is the mass per unit span, Sα = mxα b is the first moment of inertia about the elastic axis, Iα = mrα2 b2 is the moment of inertia about the elastic axis, Kh = mωh2 is the bending stiffness. xα is the airfoil static unbalance, rα is the gyration radius of the airfoil around the elastic axis, ωh is the uncoupled natural frequency in the plunging DOF, ρ is the air density, b is the half-chord length. As usual, the elastic axis and aerodynamic center of the airfoil are located at the distances of ab and xac b after the mid-chord point, respectively. In Eq. (1), M (α ) represents the nonlinear structural restoring moment. When M (α ) = Kα α , the airfoil is degenerated to a linear structural model, where Kα = Iα ωα2 is the torsion stiffness and ωα is the uncoupled natural frequency in the pitching DOF. In the present study, the free-play nonlinearity is considered, which can be expressed as

⎧ Kα (α − δ ) α ≥ δ ⎪ M (α ) = ⎨ 0 −δ<α<δ ⎪ ⎩ Kα (α + δ ) α ≤ − δ

(2)

where δ denotes the measurement of free-play. Introducing the mass ratio, μ = m /πρb2, Eq. (1) can be written as

⎧ ⎛ ¨ ⎞ ⎛ ⎞2 2 ⎪ 1 ⎜ h + xα α¨ ⎟ + ⎜ ωh ⎟ h = U ( − cl ) ⎜ ⎟ ⎝ ⎠ 2 ⎪ ωα ⎝ b ωα b πμ ⎠ ⎨ ⎞ ⎪ 1 ⎛ h¨ M (α ) U2 ⎜ xα + rα2 α¨ ⎟ + rα2 = (2cm ) ⎪ 2 ⎪ πμ b K ω ⎝ ⎠ α ⎩

(3)

α

where U = V /bωα is the speed index. For this aeroelastic system, {h/b displacements and the generalized aerodynamic forces, respectively.

α}T

and { − cl 2cm

}T

are served as the generalized

2.2. Aerodynamic nonlinearity modeling Based on the authors' previous work (He et al., 2014), the aerodynamic describing function is used to model the nonlinear aerodynamics in transonic air flow in the present study. And Euler equations are adopted to calculate the nonlinear aerodynamic forces. In the aerodynamic describing function method, Euler equations are regarded as an implicit nonlinear

4

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

aerodynamic system with the generalized displacements as inputs and the generalized aerodynamic forces as outputs. And the describing function is the ratio of the first order harmonic component of the complex output and the harmonic input. According to the definition of describing function, the input of CFD solver, i.e. pitching (or plunging) motion of the airfoil, is taken as a sinusoid function,

α (t ) = α0 sin (ωt )

(4)

where α0 denotes the amplitude of LCO and ω is the frequency. And the corresponding output of the CFD solver, aerodynamic moment coefficient, for instance, is periodic which can be expanded by Fourier series. 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 coefficients of (cm )c and (cm )s can be obtained by fitting the time history of aerodynamic moment coefficient to the above equation. And the aerodynamic describing function can be expressed as

Q mα =

(cm ) cs iφ (cm ) cs e = [cos (φ) + isin (φ)] = R (Q mα ) + iI (Q mα ) α0 α0

(6)

The detailed process to build the aerodynamic describing function can be found in the work of He et al. (2014), and the generalized aerodynamic forces corresponding to the airfoil motion can be written as

⎤⎧ ⎫ ⎧ ⎫ ⎡ ⎧ h⎫ ⎪ − cl ⎪ ⎢ Q lh Q lα ⎥ ⎪ h ⎪ ⎪ ⎪ ⎨ ⎨ b ⎬ = [Q ] ⎨ b ⎬ ⎬= 2cm ⎪ ⎢ Q mh Q mα ⎥ ⎪ ⎪ ⎪ ⎪ ⎦⎩ α ⎭ ⎩ ⎭ ⎣ ⎩ α⎪ ⎭

(7)

where [Q ] is the aerodynamic describing function, which is related to Mach number, reduced frequency k = ωb /V and the amplitude of LCO. When building the aerodynamic describing function, if the amplitudes of motion, namely α0 and h0 /b, are set to be sufficiently small in CFD calculations, the aerodynamic describing functions are reduced to linear aerodynamic influence coefficients, which are only related to Mach number and reduced frequency. Thus, we can obtain the dynamically linear aerodynamics, and only statically nonlinear aerodynamics is taken into account. When the amplitudes of motion are set to be specific LCO amplitudes, we can obtain the corresponding nonlinear aerodynamics. 2.3. Structural nonlinearity modeling In our previous study (Yang and Zhao, 1988), the harmonic balance analysis, which can also be called structural describing function method, was used to examine the LCO behavior of an aeroelastic airfoil with nonlinear pitching stiffness. The structural describing function method is a computationally efficient approach for modeling nonlinear dynamic systems when the response is periodic in time. A fundamental harmonic solution, i.e. Eq. (4), for the pitching motion is firstly assumed. Then from Eq. (2) and the structural describing function approach, the linearized equivalent stiffness can be expressed as

Keq = feq Kα

(8)

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

feq

⎧0 α0 < δ ⎪ ⎪ ⎫ ⎧ ⎛ ⎞ ⎛ δ ⎞2 ⎪ =⎨ 2⎪ δ δ 1 − ⎜ ⎟ ⎬ α0 ≥ δ ⎪ 1 − ⎨ arcsin ⎜⎜ ⎟⎟ + ⎝ α0 ⎠ ⎪ π⎩ ⎪ ⎪ ⎝ α0 ⎠ α0 ⎭ ⎩

(9)

So the structural restoring moment corresponding to the oscillation of airfoil can be obtained,

M (α ) = feq Kα α

(10)

The structural describing function method has some limitations as well. Gelb and Vander Velde (1968) noted that the fundamental limitation is that the form of the input signal of the nonlinearity must be guessed in advance and the analysis only answers specific questions asked by the assumed input. As shown in Eq. (4), a sinusoid function is assumed as the input of the nonlinearity, so it is not surprising that the magnitude of the LCO motion that is predominantly period-1 for an aeroelastic airfoil can be predicted well by using the structural describing function reported by Price et al. (1995). It is believed that the describing function method at least permits the qualitative LCO behavior of the nonlinear aeroelastic system.

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

5

2.4. Flutter solution When the airfoil is vibrating harmonically, this is always the case in critical flutter condition,

⎧ h ⎫ ⎧ h0 ⎫ ⎪ ⎪ ⎪ ⎪ ⎨ b ⎬ = ⎨ b ⎬ eiωt ⎪ ⎭ ⎩ α0 ⎪ ⎭ ⎪ ⎩ α⎪

(11)

Substituting Eqs. (7), (10) and (11) into Eq. (3) yields the flutter equation in frequency domain

⎤⎧ ⎫ ⎧ ⎫ ⎧ ⎫ ⎡ ⎛ ⎞2 ω ⎪ 0 ⎥⎪ ⎪ h⎪ ⎪ U2 ⎪ h⎪ ⎪ ⎪ h⎪ ⎪ ⎢⎜ h⎟ ⎛ ω ⎞2 ⎡ 1 xα ⎤ ⎪ ⎥⎨ b ⎬ = ⎥ ⎨ ⎬ + ⎢ ⎝ ωα ⎠ −⎜ ⎟ ⎢ [Q ] ⎨ b ⎬ ⎝ ωα ⎠ ⎣ xα rα2 ⎦ ⎪ b ⎪ ⎢ πμ ⎥⎪ α ⎪ ⎪ α⎪ α ⎪ ⎪ rα2 feq ⎥⎦ ⎪ ⎩ ⎪ ⎭ ⎩ ⎪ ⎭ ⎩ ⎪ ⎭ ⎢⎣ 0

(12)

where ω/ωα is called the frequency ratio. Adding the artificial structural damping into Eq. (12), the V-g method (Hodges and Pierce, 2011) is used to obtain the flutter speed index and flutter frequency ratio. It is important to note that if the shock motion is small and in proportion to the structural motion, the dynamically linear aerodynamics can be used to model the unsteady aerodynamic forces according to Dowell et al. (2004b). Moreover, by setting feq = 1, the above flutter solution is also suitable for flutter analysis with linear structural model. With the dynamically linear aerodynamics and linear structural model, we can obtain the linear critical flutter speed index Ul and the linear flutter frequency ratio ωl /ωα . 2.5. Stability of limit cycle oscillation The study of the stability of LCO is of great importance because the unstable LCO can not 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.

3. The aeroelastic model The aeroelastic airfoil studied here is taken from the work of Kousen and Bendiksen (1994), in which a free-play nonlinearity is assumed in the pitching DOF. The airfoil of this model is NACA 64A010, and the relevant parameters are xα = 0.2, rα2 = 0.29, μ = 60, ωh/ωα = 0.34335, a = − 0.2, δ = 1°. The computed Mach number is 0.87. 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. 3. 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 for a NACA 64A010 airfoil with a sinusoid oscillating motion for α0 = 1° and k ¼0.1 at Mach 0.8 are conducted. The time histories of the unsteady aerodynamic coefficients with different meshes are shown in Fig. 4. It can be seen that the amplitude of cl obtained by using Mesh 1 is a little smaller than those obtained by using Mesh 2 and Mesh 3, respectively, and the amplitude of cm obtained by using Mesh 1 is a little larger than those obtained by using Mesh 2 and

Fig. 3. Computational grids for NACA 64A010 airfoil: (a) Overall. (b) Close-up.

6

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

Mesh 1 Mesh 2 Mesh 3

0.2

0.03

0.15

Mesh 1 Mesh 2 Mesh 3

0.02

0.1 0.01

0

cm

cl

0.05

-0.05

0

-0.01

-0.1 -0.02

-0.15 -0.2

0

0.1

0.2

0.3

Time (s)

0.4

0.5

0.6

-0.03

0

0.1

0.2

0.3

0.4

0.5

0.6

Time (s)

Fig. 4. Comparison of aerodynamic coefficients of NACA 64A010 airfoil with different meshes (Mach 0.8, α0 = 1° , k ¼ 0.1 and a = − 0.6 ): (a) Lift coefficient. (b) Aerodynamic moment coefficient.

Mesh 3, respectively, but results obtained by using Mesh 2 and Mesh 3 are nearly the same. It reveals that both Mesh 2 and Mesh 3 can be used to calculate the unsteady aerodynamic forces. Considering the computational costs, Mesh 2 is adopted in this section. A brief introduction of the procedure to run a nonlinear flutter analysis at a prescribed Mach number is described as follows. 1. Select a LCO amplitude. 2. Use Eq. (9) to calculate the structural describing function. 3. For several reduced frequencies, use a sinusoidal motion with the selected LCO amplitude as the input of the CFD solver, perform CFD simulation and get the first-order harmonic components of the computed generalized aerodynamic forces by curve fitting, and build the aerodynamic describing functions by using Eq. (6). 4. Solve flutter equation, i.e. Eq. (12), by using V-g method to obtain U and ω/ωα . Note that how to choose the reduced frequency in V-g method should generally undergo a process of trial and error, and in the present study we specified a set of trial reduced frequency values from 0.1 to 0.4 with a step of 0.05. Repeating the above procedures for different LCO amplitudes and depicting the curve of U versus LCO amplitude, the bifurcation diagram can be obtained as shown in Fig. 5. This figure also shows the bifurcation diagram obtained by using the time marching approach based on Euler equations (Kousen and Bendiksen, 1994). It should be mentioned that the reduced velocity in Kousen and Bendiksen (1994) is multiplied by 2 to obtain speed index U defined in the present study for comparisons. As shown in Fig. 5, when the LCO amplitude is less than 8°, the bifurcation diagram from the present method agrees well with that obtained by time marching approach (Kousen and Bendiksen, 1994). However, as the LCO amplitude increases further, differences of the results between these two methods become larger. Some explanations regarding the differences of the flow speed shown in Fig. 5 considering both aerodynamic and structural nonlinearities should be addressed. Firstly, the linear flutter speed index is a little different, nearly 4 from Kousen and Bendiksen (1994) and 3.75 from the present method. In conventional numerical calculations only considering free-play nonlinearity, the LCO amplitude will be very large or nearly exhibit divergent oscillations when flow speed approaches the linear flutter speed. That is to say, the asymptote of the curve for LCO amplitude versus flow speed should be vertical to the speed axis at the linear flutter speed. So the difference of the linear flutter speed can cause bifurcation diagram moving along the speed axis. Secondly, the accuracy of the aerodynamic describing function method decreases as the LCO amplitude increases, which is caused by the curve fitting of aerodynamic force coefficients. It can be illustrated by Fig. 6 that the difference between the curve fitting and direct CFD simulation results becomes evident as the LCO amplitude increases. Last but not least, the accuracy of superposition principle used in building the aerodynamic describing function also decreases as the LCO amplitude increases (He et al., 2014). However, the LCO behavior obtained by using the present method in general is in reasonable agreement with that obtained by using the time marching approach. Considering that we pay more attention on the LCO behaviors of an

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

7

Kousen and Bendiksen (1994) Nonliear aerodynamics + nonlinear structure Nonlineare aerodynamics + linear strucuture Dynamically linear aerodynamics + nonliear structure

14

LCO amplitude (°)

12 10 8 6 4 2 0 0.5

1

1.5

2

2.5

3

3.5

4

4.5

U Fig. 5. LCO amplitude versus speed index for NACA 64A010 airfoil at Mach number 0.87.

aeroelastic airfoil in the present study, it is feasible to use the aerodynamic describing function method for transonic LCO solutions. To evaluate the contribution of aerodynamic nonlinearity to the LCO behaviors, the results of the linear structural model with nonlinear aerodynamics are also illustrated in Fig. 5. It is observed that the curve of LCO amplitude versus airspeed is nearly perpendicular to the airspeed axis, indicating a so-called linear LCO behavior trend as mentioned in Section 1. It reveals that the aerodynamic nonlinearity is weak at Mach 0.87. Therefore, it is not surprising to observe that the LCOs obtained by dynamically linear aerodynamics and nonlinear aerodynamics with nonlinear structural model are in good agreements, as shown in Fig. 5. From the above results, it is found that the assumption of dynamically linear aerodynamics is suitable to obtain the LCO behavior at this Mach number. The dynamically linear aerodynamics will be assumed in the following sections to explore the LCO behaviors at different Mach numbers. Note that this type of aerodynamics was also applied in transonic LCO analysis for an aeroelastic airfoil with a structural nonlinearity by Dowell et al. (2004b) and Jones et al. (2007). It is important to note that though the aerodynamics is assumed to be dynamically linear, the structural model is nonlinear in the present study. So the limiting mechanism for the LCO in the following sections is actually caused by freeplay structural nonlinearity.

4. Results and discussions 4.1. LCO behaviors By using the structural describing functions and dynamically linear aerodynamics, the bifurcation diagrams are obtained from Mach 0.3 to 0.8 with a step increment of 0.1 and from Mach 0.81 to 0.95 with a step increment of 0.01. Representative bifurcation diagrams at several different Mach numbers are shown in Fig. 7. In the Mach number range from 0.3 to 0.83, a type of subcritical bifurcation is observed, which consists of both stable branch for large amplitudes and unstable branch for small amplitudes. When the Mach number ranges from 0.84 to 0.91, the unstable branch turns to be stable. And then the stable branch for small LCO amplitudes turns back to be unstable again at Mach numbers ranging from 0.92 to 0.95. It is clear that the LCO behaves differently in different Mach number ranges, and the major difference is that the LCOs can be achieved even at very small flutter speeds in the Mach number range from 0.84 to 0.91. Several LCO results (predominantly period-1) from time marching simulation based on Euler equations are also plotted in Fig. 7. There are reasonably good agreements between results of the present method and the time marching approach, demonstrating the validity of our method for qualitative LCO behavior 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. 7 as well. Redrawing the LCO results of Figs. 7 and 8 shows the flutter boundary for different LCO amplitudes. It can be seen from

8

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

CFD Curve fitting

CFD Curve fitting

0.08

0.015

0.06

0.01

0.04

0.005

0.02

Cm

Cm

0.02

0

0

-0.005

-0.02

-0.01

-0.04

-0.015

-0.06

-0.02

-0.08 0.4

0.42

0.44

0.46

0.48

0.5

0.52

0.54

0.4

0.42

0.44

time (s)

0.46

0.48

0.5

0.52

0.54

time (s)

CFD Curve fitting

0.1 0.08 0.06

Cm

0.04 0.02 0 -0.02 -0.04 -0.06 -0.08 -0.1

0.4

0.42

0.44

0.46

0.48

0.5

0.52

0.54

time (s)

Fig. 6. Comparison of the aerodynamic moment coefficients obtained by CFD simulation and curve fitting at different LCO amplitudes (k¼ 0.2, Mach 0.87): (a) α0 = 2°, (b) α0 = 8° , (c) α0 = 10° .

Fig. 8(a) that, in the narrow Mach number range from 0.84 to 0.91, the LCOs with different amplitudes are observed over a wide range of flutter speed. And this map of LCO regions looks like the “chimney” phenomenon reported by Silva et al. (2000) and Edwards et al. (2001). From Fig. 8(a), the LCO amplitudes inside the chimneys, varied from 1° to 1.2°, are insensitive to the flutter speed, which is consistent with the conclusion on the chimney phenomenon from Bendiksen (2015). Furthermore, as shown in Fig. 8(b), the flutter frequencies change greatly with the variation of Mach number from 0.8 to 0.95 indicating a flutter mode change, which is also a remarkable feature of the chimney phenomenon according to Dowell et al. (2004b). Returning to our topic of LCO behaviors for different Mach numbers, we will address why the unstable LCOs with small amplitudes turn to be stable ones in the narrow Mach number range and how to estimate the low and high Mach numbers of this range in the following sections. Comparing Fig. 7 with Fig. 8, we can see that this specific Mach number range will be the same as that of the chimney region, so the low and high Mach numbers of this region will also be called as the left and right boundaries, respectively, of the chimney for convenience. Furthermore, the reasons for the LCO trends behaving variously at different Mach number ranges would be obviously the same as those for the chimney phenomenon in the present study. Thus, we could get started from the explanations for the chimney phenomenon.

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

Present, Ma 0.8 Present, Ma 0.83 Present, Ma 0.84 Present, Ma 0.87 Present, Ma 0.91 Present, Ma 0.92 Present, Ma 0.95 Time-marching, Ma 0.8 Time-marching, Ma 0.87 Time-marching, Ma 0.95

5

4

LCO amplitude (°)

9

3

2 Stable LCO

1 Unstable LCO Stable LCO

0

0

1

2

3

4

U Fig. 7. LCO amplitude versus speed index for NACA 64A010 airfoil at different Mach numbers. “Time marching” stands for time marching approach based on Euler equations.

Fig. 8. Speed index and frequency ratio verses Mach number for NACA 64A010 airfoil with different LCO amplitudes (the LCO amplitude is marked on the line).

4.2. Mach number freeze and the high Mach number prediction According to Bendiksen (2015), Mach number freeze phenomenon or transonic stabilization law plays an important role in determining the shape of the chimney regions. Originally, the Mach number freeze is a phenomenon that occurs near Mach 1 in two-dimensional steady flow, where the local Mach number at the points on or near the airfoil surface ahead of the shock “freezes” and becomes independent of the free stream Mach number. And a corresponding unsteady stabilization law for transonic air flow was proposed by Bendiksen (2012), which leads to a stabilization law for flutter. According to Bendiksen (2012), the Mach number freeze begins when the shocks on the upper and lower surface have reached the trailing edge of the airfoil in two dimensional air flow, and the specific Mach number is called the freeze Mach number.

10

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18 -1

Ma = 0.8: 0.95 with step of 0.01

Mach 1.3 1.2 1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

Pressure coefficient

-0.5

0

Ma = 0.91

0.5

1

1.5

0

0.2

0.4

x/c

0.6

0.8

1

Mach

Mach

1.3 1.2 1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

1.3 1.2 1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

Fig. 9. Steady flow field for NACA 64A010 airfoil at different Mach numbers at zero angle of attack: (a) Steady pressure coefficient distribution. (b) Mach number contours at Mach 0.91. (c) Mach number contours at Mach 0.92. (d) Mach number contours at Mach 0.93.

For the present NACA 64A010 airfoil, Fig. 9 shows the steady pressure coefficient distributions and Mach number contours at different Mach numbers at zero angle of attack. And it is evident that the shock forms and moves towards the trailing edge as the Mach number increases. The shock wave has reached the trailing edge when the Mach number is increased to 0.92, indicating the freeze Mach number. Bendiksen (2015) employs the Mach number freeze phenomenon to explain the chimney phenomenon as follows. The amplitude and phase of the unsteady aerodynamic forces reverse immediately before the onset of the freeze Mach number. These reversals of the aerodynamic forces will trigger the LCOs dominated by the single degree of freedom (SDOF) flutter, and a strong quenching of the LCOs follows at the freeze Mach number. When this phenomenon occurs, weak unstable LCO can persist down to the low dynamic pressure, thus creating a deep chimney in the transonic flutter boundary. The Mach number freeze phenomenon is also adopted to explore our observations herein. As shown in Fig. 7, we can infer that as the Mach number is increased to 0.84, i.e. the left boundary (low Mach number) of the chimney region, the amplitude and phase reversals of the unsteady aerodynamic forces will happen and trigger the SDOF flutter. Furthermore, these reversals turn the unstable LCOs with small amplitudes in subsonic regime to stable ones at very low flutter speeds. When the Mach number is increased to 0.92, i.e., the freeze Mach number, the quenching of LCOs caused by the freeze phenomenon occurs, which forms the right boundary (high Mach number) of the chimney. Meanwhile, the stable LCO branch for small amplitude turns back to unstable and the flutter speed increases abruptly to a very high speed, as shown in Fig. 8. So the Mach number freeze phenomenon can provide a physics-based explanation for the chimney phenomenon as well as for that the LCOs behave variously at different Mach number ranges. The high Mach number can be well predicted by the freeze Mach number, but how to estimate the low Mach number remains unsolved, which will be addressed in the following sections.

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

11

4.3. Single degree of freedom flutter Since the SDOF flutter is of great importance in the chimney phenomenon, a detailed discussion on the SDOF flutter for the present aeroelastic system will be conducted in this section. Only the pitching motion is retained in Eq. (12), and the motion equation of SDOF flutter system in frequency domain can be written as

⎛ ω ⎞2 2U 2 Q mα α −⎜ ⎟ α + feq α = ⎝ ωα ⎠ πμrα2

(13)

According to the procedure of V-g method (Hodges and Pierce, 2011), introducing the artificial structural damping coefficient g, the SDOF flutter solution can be obtained:

g=

feq I (Q mα ) R (Q mα ) + πμrα2 k2/2

ω/ωα =

(14a)

feq 2R (Q mα )/πμrα2 k2 + 1

(14b)

ω/ωα U= k

(14c)

For a specific reduced frequency, the system is stable if the structural damping coefficient is negative, while the system is unstable when the structural damping coefficient is positive. When the damping coefficient approaches zero, the critical reduced frequency and the flutter speed index are achieved. From Eqs. (4) to (6) and Eq. (14a), it is apparent that the phase reversal of the unsteady aerodynamic force can reflect the transition of damping coefficient from negative to positive. It is important to note that the phase reversal of the unsteady aerodynamic force, rather than both of the amplitude and phase reversals according to Bendiksen (2015), is equivalent to the emergence of SDOF flutter for the present aeroelastic system. 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. It means that at a specific Mach number, the critical reduced frequencies of the current SDOF system are all the same for each prescribed LCO amplitude. Thus, we take the aeroelastic system with linear structural model as an example to illustrate the specific Mach number at which the SDOF flutter may occur. Fig. 10 shows the curves of structural damping coefficient versus speed index between Mach 0.81 to Mach 0.95. And it can be found that the SDOF flutter may occur only in the Mach number range from 0.84 to 0.91, which agrees well with the Mach number range of the chimney region. It should be mentioned that the following important traits are all shared by the flutter characteristics in the chimney region, but we only take the result at Mach 0.87 as a representative example in the following discussion. To confirm the emergence of the SDOF flutter occurred at the chimney region, the flutter speeds of the SDOF system only considering pitching DOF and the original two DOF aeroelastic system considering both pitching and plunging DOFs are compared, as shown in Fig. 11. From this figure, we can see that when the LCO amplitude is less than 1.1°, the flutter speeds 0.6 Ma = 0.81 Ma = 0.83 Ma = 0.84 Ma = 0.87 Ma = 0.91 Ma = 0.92 Ma = 0.95

0.4

g

0.2

0

-0.2

-0.4

-0.6

2.5

3

3.5

4

4.5

5

5.5

6

U Fig. 10. Speed index versus structural damping coefficient with linear structure for different Mach numbers.

12

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

4

LCO amplitude (°)

3

Flutter speed for 2 DOF aeroelstic system 2

Flutter speed for SDOF aeroelstic system α0 = 1.1°

1

0

0

1

2

3

4

U Fig. 11. LCO amplitude versus speed index for single degree of freedom and two degree of freedom aeroelastic systems for different LCO amplitudes at Mach 0.87.

of the SDOF system are almost the same as those of the two DOFs aeroelastic system. Such results indicate that the LCO with small amplitude at low flutter speed is dominated by the SDOF flutter. To identify the flutter type of the original two DOFs aeroelastic airfoil, the flutter frequencies and the natural frequencies for different LCO amplitudes are compared at Mach numbers of 0.5, 0.87 and 0.95, as shown in Fig. 12. For the Mach number within the chimney region, for example, at Mach 0.87, the flutter frequency corresponding to the amplitude ranging from 1.0° to 1.1° is nearly the same as the first natural frequency of the system, so the flutter type of aeroelastic airfoil with the small amplitude should be SDOF flutter according to Bendiksen (2015). When the LCO amplitude increases, the flutter frequency becomes apparently larger than the first natural frequency and lies between the first and second natural frequencies of the system, indicating that the effect of the second natural mode on the flutter speed becomes significant. Such result reveals that the flutter with the large LCO amplitude can be regarded as the coupled bending-torsion flutter. For the Mach number outside of the chimney, for instance, at Mach 0.5 and Mach 0.95, the flutter mode should always be the 1

ω2/ωα 0.8

Ma = 0.87

0.6

ω/ωα

Ma = 0.95 Ma = 0.5

0.4

ω1/ωα 0.2

0

1

1.5

2

2.5

3

LCO amplitude (°) Fig. 12. Comparisons of the frequency ratio for different LCO amplitudes at different Mach numbers (ω1/ωα and ω2/ωα are the first and second natural frequency ratios, respectively).

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

13

0.008

0.18 Ma = 0.81 Ma = 0.87 Ma = 0.95

0.16

Ma = 0.81 Ma = 0.87 Ma = 0.95

0.006

0.14

0.004

0.12 0.002

cl

cm

0.1

0

0.08 0.06

-0.002

0.04 -0.004

0.02 0 0.1

-0.006 0.2

0.3

Angle of attack (°)

0.4

0.5

0.1

0.2

0.3

0.4

0.5

Angle of attack (°)

Fig. 13. Lift coefficient and moment coefficient versus angle of attack for NACA 64A010 airfoil at different Mach numbers: (a) Lift coefficient. (b) Aerodynamic moment coefficient.

coupled bending-torsion mode no matter whether the LCO is stable or not, which can also be observed from Fig. 12. From the above discussion, we can see that when the Mach number is increased to the chimney region, the phase reversal of the unsteady aerodynamic forces will trigger the SDOF flutter in pitching motion. Meanwhile, the unstable LCO with small amplitude, caused by the coupled bending-torsion mode, turns to be a stable one dominated by the SDOF flutter. 4.4. The low Mach number prediction As mentioned in Section 4.3, the flutter type within the chimney region is the SDOF flutter, so a possible answer to the question of how to estimate the low Mach number of the chimney region may be found from the condition for the occurrence of the SDOF flutter. Smilg (1949) noted that the SDOF flutter occurs in subsonic incompressible flow if the airfoil moment of inertia is sufficiently large and the elastic axis of the airfoil is located at a point that is ahead of the airfoil quarter-chord but not too far ahead of the airfoil leading edge. According to the thin airfoil aerodynamic theory, the aerodynamic center is located at the quarter-chord point in the subsonic incompressible flow. It reminds us whether the conditions for the occurrence of the SDOF flutter is related to locations of the elastic axis and aerodynamic center. Fig. 13 shows the lift coefficient and moment coefficient curves at low angles of attack for different Mach numbers. It is evident that the lift coefficient and moment coefficient vary linearly with the angle of attack of the airfoil, indicating that the aerodynamic center exists as a fixed point on the airfoil according to Anderson (2007). And the location of the aerodynamic center can be determined by the following equation:

x ac = − 2

dcm dcl / +a dα dα

(15)

As in Smilg (1949), the location of the elastic axis (a) and the reduced velocity (V /ωb or 1/k ) required to satisfy g ¼0 are plotted in one figure, then the condition for the occurrence of the SDOF flutter is obtained. Similarly, the curve of a versus V /ωb for NACA 64A010 airfoil at three typical Mach numbers are plotted in Fig. 14. It can be seen that there should be a certain relationship between the conditions for the occurrence of the SDOF flutter and locations of the elastic axis and aerodynamic center. Within scope of the present study, this relationship can be described as that the elastic axis should be located near or even ahead of the aerodynamic center, then the SDOF flutter is possible to happen. It should be noted that the effect of the airfoil moment of inertia on SDOF flutter is not discussed. Generally speaking, with the increasing of Mach number, the aerodynamic center of the airfoil will move from the quarter-chord point (xac = − 0.5) in subsonic air flow to the middle-chord point (xac = 0.0) in supersonic air flow according to the basic concept of the thin airfoil aerodynamic theory. For the present aeroelastic airfoil (a = − 0.2), the elastic axis is located ahead of the aerodynamic center in the subsonic air flow, where the SDOF flutter seems to be impossible to occur. But as Mach number increases, the aerodynamic center moves towards the trailing edge, and when it lies behind the elastic axis, the SDOF flutter may occur. Therefore, it is speculated that the Mach number, at which the aerodynamic center lies near the elastic axis, can be used to estimate the low Mach number of the chimney region herein. Fig. 15 shows the locations of aerodynamic center varying with the Mach numbers. It is shown that in subsonic and high

14

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

Fig. 14. Region of instability in single degree of freedom aeroelastic system at different Mach number: (a) Mach 0.7. (b) Mach 0.83. (c) Mach 0.87.

subsonic air flow (Ma < 0.7), the aerodynamic center is nearly fixed at the quarter-chord point (xac = − 0.47). While in transonic regime, as the Mach number increases, the location of the aerodynamic center firstly moves towards the trailing edge and then reversely moves towards the leading edge. When the Mach number is increased to the freeze Mach number, i.e. 0.92, the aerodynamic center will be nearly fixed at xac = − 0.13 as shown in Fig. 15. It is important to note that if the Mach number is greater than 0.83, the aerodynamic center lies behind the elastic axis for the present aeroelastic airfoil. At this condition the SDOF flutter may occur, which implies that the low Mach number of the chimney region would be about 0.84. Such result agrees well with the LCO analysis in Section 4.1, therefore, the low Mach number of the chimney region can be well predicted by the specific Mach number at which the aerodynamic center is located near the elastic axis. From the above discussion, it is found that the characteristics of aerodynamic center movement with the increasing of Mach number plays an important role on the appearance of the SDOF flutter in the two-dimensional aeroelastic system. In the three-dimensional wing model, Bendiksen (2015) noted that LCOs dominated by a natural mode, which is also called natural mode instability, occur in the chimney region. Natural mode instability can be dominated by the first torsional mode or second bending mode for different types of wing configurations according to Bendiksen (2015). One may imagine that LCO dominated by the first torsional mode in three-dimensional wing model is similar to the SDOF flutter of pitching DOF in an aeroelastic airfoil. Then it can be speculated that the relative location between aerodynamic center of the wing and the

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

15

1

( -0.13, 0.92)

0.9

(0.26, 0.89) ( -0.2, 0.833)

0.8

Ma

0.7 0.6 0.5

(-0.47, 0.4)

0.4

Elastic axis

0.3

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

xac

0.4

0.6

0.8

1

Fig. 15. Aerodynamic center locations at different Mach numbers for NACA 64A010 airfoil.

nodal line of the first torsional mode may be related to the appearance of LCOs dominated by the first torsional mode. Keeping that thought in mind, the Mach number which causes the aerodynamic center to align with the nodal line of the torsional mode might be used to estimate the lower Mach number of the chimney region. Moreover, the analysis of the effect of wing configuration on the aerodynamic center movement with the increasing of Mach number may be helpful to determine what types of wing configurations would be susceptible to natural mode instabilities dominated by the first torsional mode. 4.5. Influence of angle of attack and viscous effects In essence, dynamically linear aerodynamics modeling is a linearization problem at the CFD steady (equilibrium) solution that can be affected by the angle of attack and viscous effects. In this section, we will explore how sensitive the bifurcation behavior is to the equilibrium point. Fig. 16 shows the typical bifurcation diagrams at different angles of attack, 0.1° and 0.5°, in inviscid air flow. Compared with Fig. 7, it is obvious that at the same Mach number, the bifurcation diagrams at angles of attack (0.1° and 0.5°) are 5

5 Ma 0.8 Ma 0.83 Ma 0.84 Ma 0.87 Ma 0.91 Ma 0.92 Ma 0.95

4

LCO amplitude (°)

LCO amplitude (°)

4

3

2 Stable LCO

1

3

2 Stable LCO

1 Unstable LCO

Stable LCO

0

Ma 0.8 Ma 0.83 Ma 0.84 Ma 0.87 Ma 0.91 Ma 0.92 Ma 0.95

0

1

2

U

Unstable LCO

Stable LCO

3

4

0

0

1

2

3

U

Fig. 16. LCO amplitude versus speed index for NACA 64A010 airfoil at different Mach numbers at the angle of attack: (a) 0.1°. (b) 0.5°.

4

16

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

5 Ma 0.8 Ma 0.83 Ma 0.84 Ma 0.87 Ma 0.88 Ma 0.89 Ma 0.95

LCO amplitude (°)

4

3

2 Stable LCO

1

Unstable LCO

Stable LCO

0

0

1

2

3

4

U Fig. 17. LCO amplitude versus speed index for NACA 64A010 airfoil at different Mach numbers considering viscous effects.

almost the same as that at zero angle of attack, except that the linear flutter speed is a little different. As we know the freeze Mach number and the Mach number at which the aerodynamic center of the airfoil lies near its elastic axis do not change significantly at such small angles of attack. Thus, the bifurcation diagram, the low and high Mach numbers of the chimney region are unchanging at small angle of attack except the linear flutter speed. Fig. 17 shows the representative bifurcation diagrams in viscous air flow. When building the dynamically linear aerodynamics, the shear stress transport (SST) turbulence model is adopted in the CFD solver. The linear flutter speed considering viscous effects is a little different from that in inviscid air flow, which is consistent with the conclusion of Prananta et al. (1998). Compared with Fig. 7, the low Mach number of the chimney region remains 0.84 but the high one changes to 0.88 when taking into account viscous effects. The LCO behaviors between Mach 0.89 and 0.91 are different from those in inviscid air flow. From the steady pressure coefficient distribution shown in Fig. 18, it is found that the freeze Mach number is still around Mach 0.92, which can also give a reasonably good prediction for the high Mach number of the chimney region in the viscous air flow. -1

Ma = 0.8: 0.95 with step of 0.01

Pressure coefficient

-0.5

0

Ma = 0.91

0.5

1

1.5

0

0.2

0.4

0.6

0.8

1

x/c Fig. 18. Steady pressure coefficient for NACA 64A010 airfoil at different Mach numbers at zero angle of attack considering viscous effects.

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

17

5. Conclusions Euler equations are employed to calculate the unsteady aerodynamic forces for transonic LCO analysis, and an aeroelastic airfoil with free-play in pitching DOF is adopted in the current study. Aerodynamic describing function and structural describing function are used to deal with aerodynamic and structural nonlinearities, respectively. Then the flutter speed index and flutter frequency ratio are obtained by V-g method. Without considering the free-play nonlinearity and only considering the aerodynamic nonlinearity, a so-called linear LCO behavior is observed. For the nonlinear structural model, the LCO solutions obtained with dynamically linear aerodynamics agree well with those obtained with nonlinear aerodynamics, indicating that the effect of aerodynamic nonlinearity is weak for this model. The dynamically linear aerodynamics is assumed to model the transonic aerodynamics of the aeroelastic airfoil for LCO analysis at different Mach numbers. When the Mach number is less than 0.83, the phenomenon of subcritical bifurcation is observed, which consists of both stable and unstable branches. However, if the Mach number increases from 0.84 to 0.91, the unstable branch for small LCO amplitudes turns to be stable. With the further increasing of Mach number, the stable branch turns back to be unstable. Specially, the LCOs can be achieved even at very small flutter speeds at Mach numbers ranging from 0.84 to 0.91. Furthermore, the flutter boundary shows that the LCOs with different amplitudes are observed over a wide range of flutter speed at this Mach number range, and the so-called chimney phenomenon is observed. The Mach number freeze phenomenon offers a physics-based explanation for the chimney phenomenon and why the unstable LCOs with small amplitudes turns to be stable in this specific Mach number range. The freeze Mach number gives a good prediction for the high Mach number of the chimney region. As we know, with the increasing of Mach number, the aerodynamic center of the airfoil moves from the quarter-chord point in subsonic air flow to the middle-chord point in supersonic air flow. In the current study, the specific Mach number, at which the aerodynamic center is located near the elastic axis, is suggested to estimate the low Mach number of the chimney region. For the present aeroelastic airfoil, as the Mach number increases, when the aerodynamic center is located behind the elastic axis, the phase reversal of the aerodynamic forces takes place and triggers the single degree of freedom flutter at a low flutter speed with small LCO amplitude in the pitching DOF. Meanwhile, the unstable LCO with small amplitude, caused by the coupled bending and torsion modes, turns to be a stable one dominated by the single degree of freedom flutter. Then if the freeze Mach number is reached, the quenching of LCOs occurs, which implies that the high Mach number of the chimney region is achieved.

Acknowledgments This work is supported by the National Natural Science Foundation of China (Grant no. 11472216 and 11102162), and the 111 project of China (Grant no. B07050). The authors wish to acknowledge Prof. Earl H. Dowell in Duke University for his comments and suggestions that improved this work.

References Anderson Jr, J.D., 2007. Fundamentals of Aerodynamics, 4th edition. The McGraw-Hill companies, New York, USA. Bendiksen, O.O., 2008. Transonic limit cycle flutter of high-aspect-ratio swept wings. J. Aircr. 45, 1522–1533. Bendiksen, O.O., 2012. Transonic stabilization laws for unsteady aerodynamics and flutter. In: 53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference No. AIAA-2012-1718, American Institute of Aeronautics and Astronautics, Honolulu, Hawaii, April. Bendiksen, O.O., 2015. Transonic flutter characteristics of advanced fighter wings. In: 56th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference No. AIAA-2015-0438, American Institude of Aeronautics and Astronautics, Kissimmee, Florida, January. Chen, T., Xu, M., Xie, L., 2014. Aeroelastic modeling using geometrically nonlinear solid-shell elements. AIAA J. 52, 1980–1993. Dowell, E.H., Clark, R., Cox, D., Curtiss, J.H.C., Edwards, J.W., Hall, K.C., Peters, D., Scanlan, R., Simiu, E., Sisto, F., Strganac, T.W., 2004a. A Modern Course in Aeroelasticity, 4th edition. Kluwer Academic Publishers, Dordrecht. Dowell, E.H., Thomas, J.P., Hall, K.C., 2004b. Transonic limit cycle oscillation analysis using reduced order aerodynamic models. J. Fluids Struct. 19, 17–27. Edwards, J.W., Schuster, D.M., Spain, C.V., Keller, D.F., Moses, R.W., April 2001. MAVRIC flutter model transonic limit cycle oscillation test. In: 42nd AIAA/ ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference. No. AIAA-2001-1291. American Institute of Aeronautics and Astronautics, Seattle, WA. 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 J. 52, 1393–1403. Hodges, D.H., Pierce, G.A., 2011. Introduction to Structural Dynamics and Aeroelasticity. Cambridge University Press, New York. Jones, D.P., Roberts, I., Gaitonde, A.L., 2007. Identification of limit cycles for piecewise nonlinear aeroelastic systems. J. Fluids Struct. 23, 1012–1028. Kholodar, D.B., Dowell, E.H., Thomas, J.P., Hall, K.C., 2004. Limit-cycle oscillation of a typical airfoil in transonic flow. J. Aircr. 44 (5), 1067–1072. Kim, D.-H., Lee, I., 2000. Transonic and low-supersonic aeroelastic analysis of a two-degree-of-freedom airfoil with a freeplay non-linearity. J. Sound Vib. 234, 859–880. Kousen, K.A., Bendiksen, O.O., 1994. Limit cycle phenomena in computational transonic aeroelasticity. J. Aircr. 31, 1257–1263. Munteanu, S., Rajadas, J., Nam, C., Chattopadhyay, A., 2005. Reduced-order-model approach for aeroelastic analysis involving aerodynamic and structural nonlinearities. AIAA J. 43, 560–571. Prananta, B., Hounjet, M., Zwaan, R., 1998. Two-dimensional transonic aeroelastic analysis using thin-layer Navier–Stokes method. J. Fluids Struct. 12 (6), 655–676. Price, S., Alighanbari, H., Lee, B., 1995. The aeroelastic response of a two-dimensional airfoil with bilinear and cubic structural nonlinearities. J. Fluids Struct.

18

Z. Yang et al. / Journal of Fluids and Structures 66 (2016) 1–18

9 (2), 175–193. Silva, W.A., Keller, D.F., Florance, J.R., Cole, S.R., Scott, R.C., 2000. Experimental steady and unsteady aerodynamic and flutter results for HSCT semispan models. In: 41st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference No. AIAA-2000-1697, American Institute of Aeronautics and Astronautics, Atlanta, GA, April. Smilg, B., 1949. The instability of pitching oscillations of an airfoil in subsonic incompressible potential flow. J. Aeronaut. Sci., 16. Thomas, J.P., Dowell, E.H., Hall, K.C., 2002. Nonlinear inviscid aerodynamic effects on transonic divergence, flutter and limit cycle oscillations. AIAA J. 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. J. Aircr. 41, 1266–1274. Yang, Z.C., Zhao, L.C., 1988. Analysis of limit cycle flutter of an airfoil in incompressible flow. J. Sound Vib. 123, 1–13.