Published: 30 September 2014

# The flutter of a two-dimensional lifting surface with cubic nonlinearities in a supersonic flow field

Hongli Zhang1
Fangqi Chen2
Xiaohua Zhang3
1Department of Mechanics, Nanjing University of Aeronautics and Astronautics, Nanjing, 210016, China
2partment of Mathematics, Nanjing University of Aeronautics and Astronautics, Nanjing, 210016, China
3Department of Information Technology, Jiangsu Institute of Commerce, Nanjing, 211168, China
Corresponding Author:
Hongli Zhang
Views 49

#### Abstract

The flutter of a two-dimensional lifting surfaces with cubic structural and aerodynamic nonlinearities in a supersonic flow field is considered. By using the maple program, the normal form of a general system of motion of a rigid airfoil with flight Mach number, nonlinear stiffness coefficient and isentropic gas coefficient is computed. The bifurcation theory is used to obtain the universal unfolding of the normal form, and all the dynamic behavior of the universal unfolding is discussed. In addition, the flutter boundary and its character are derived. This effective analysis method is applied to two numerical examples, and the influence of the structural and aerodynamic parameters on the character of flutter instability is investigated. Finally, the numerical simulations obtained by using fourth-order Runge-Kutta method will illustrate the quality of this analysis method.

## 1. Introduction

In recent years, many problems concerning the 2-degree-of-freedom (2-DOF) aeroelatstic system with cubic nonlinearity have been successfully studied using the center manifold theory and the normal form method [1-3]. There exist many potential sources of nonlinearities, which can have significant effect on an aeroelastic response. The nonlinearities of the aeroelastic system can be structural [4-5] or aerodynamic [6].

Breitbach [4] presented a detailed review of some possible structural nonlinearities and their effect on aeroelastically induced vibrations. Raghothama and Narayanan [7] investigated periodic oscillations and bifurcations of a two-dimensional airfoil in plunge and pitch motions with cubic pitch stiffness in incompressible flow. Shahrzad and Mahzoon [8] studied the limit cycle oscillation (LCO) flutter of a two-dimensional wing with nonlinear pitch stiffness. Liu, Wong and Lee [9] adopted a point transformation technique to study the nonlinear behavior of a two-dimensional aeroelastic system with freeplay models. Zhao and Zhang [10] investigated the LCO flutter and chaotic motion of a 2-DOF airfoil system with a combination of freeplay and cubic nonlinearities. The Krylov-Bogoliubov-Mitropolsky (KBM) method was used by Yang [11] to study the LCOs of an airfoil with the freeplay nonlinearity. Kim and Lee [12] analyzed the dynamic behavior of a flexible airfoil with the freeplay nonlinearity. Price et al. [13] discovered the LCOs at flow speeds below the linear flutter boundary in a system with a hard pitch spring. Abbas et al. [14] investigated the flutter and post-flutter of a two-dimensional double-wedge lifting surface. The behavior of bifurcation point and the influence of the system parameters on the stability and amplitudes of the LCO flutter were analyzed in references [15-16]. The relation between frequency and airfoil parameters was obtained by means of the harmonic balance method in reference [17].

In this paper, an effective analysis method to the problem of flutter of a two-dimensional airfoil with cubic structural and aerodynamic nonlinearities is proposed. This analysis method enables one to make a parametric study. Firstly, the normal form and universal unfolding are deduced by means of the maple program and the bifurcation theory. Secondly, super- and sub-critical Hopf bifurcations are analyzed. Finally, numerical examples and simulations are given to show the efficiency of this analysis method. Moreover, the influence of the structural and aerodynamic parameters, representations $B$, $M$, and $\gamma$, on the airfoil safety is studied.

Fig. 1Airfoil section

## 2. Nonlinear model of a two-dimensional lifting surface

Consider a two-dimensional lifting surface featuring plunging and twisting degrees of freedom, elastically constrained by a linear translational spring and nonlinear torsional spring, as shown in Fig. 1, exposed to a supersonic flow field. The equations of motion of the lifting surface are [18]:

1
$\begin{array}{l}m{h}^{\text{'}\text{'}}\left(t\right)+{S}_{\alpha }{\alpha }^{\text{'}\text{'}}\left(t\right)+{c}_{h}{h}^{\text{'}}\left(t\right)+{K}_{h}h\left(t\right)={L}_{\alpha }\left(t\right),\\ {S}_{\alpha }{h}^{\text{'}\text{'}}\left(t\right)+{I}_{\alpha }{\alpha }^{\text{'}\text{'}}\left(t\right)+{c}_{\alpha }{\alpha }^{\text{'}}\left(t\right)+{\stackrel{-}{M}}_{\alpha }={M}_{\alpha }\left(t\right),\end{array}$

where $h$ is the plunging displacement (positive downward), $\alpha$ the twist angle about the pitch axis (positive nose up), the primes denote differentiation with respect to time $t$, $m$ the structural mass per unit span, ${S}_{\alpha }$ the static unbalance about the elastic axis, ${I}_{\alpha }$ the mass moment of inertia about the elastic axis of the airfoil, ${\mathrm{c}}_{h}$ and ${\mathrm{c}}_{\mathrm{\alpha }}$ are linear plunging and pitching viscous damping coefficients, while ${K}_{h}$ is the plunging stiffness coefficient. And ${\stackrel{-}{M}}_{\alpha }$ represents the overall restoring moment that is connected to the pitch angle by:

2
${\stackrel{-}{M}}_{\alpha }={K}_{\alpha }\alpha \left(t\right)+{\delta }_{S}{\stackrel{^}{K}}_{\alpha }{\alpha }^{3}\left(t\right),$

where ${K}_{\mathrm{\alpha }}$ and ${\stackrel{^}{K}}_{\alpha }$ are the linear and nonlinear pitching stiffness coefficients, respectively. The tracer ${\delta }_{S}$ that identifies this type of nonlinearities can take the value 1 or 0 depending on whether the nonlinearity is accounted for or discarded, respectively. In addition, the expressions of the aerodynamic lift and moment are:

3
$\begin{array}{l}{L}_{\alpha }\left(t\right)=-\frac{b{U}_{\infty }{\rho }_{\infty }}{3{M}_{\infty }}\lambda \left\{12{U}_{\infty }\alpha \left(t\right)+{M}_{\infty }^{2}{U}_{\infty }\left(1+\gamma \right){\lambda }^{2}{\alpha }^{3}\left(t\right)+12\left[{h}^{\text{'}}\left(t\right)+\left(b-{x}_{EA}\right){\alpha }^{\text{'}}\left(t\right)\right]\right\},\\ {M}_{\alpha }\left(t\right)=\frac{b{U}_{\infty }{\rho }_{\infty }}{3{M}_{\infty }}\lambda \left\{12{U}_{\infty }\left(b-{x}_{EA}\right)\alpha \left(t\right)+{M}_{\infty }^{2}{U}_{\infty }\left(b-{x}_{EA}\right)\left(1+\gamma \right){\lambda }^{2}{\alpha }^{3}\left(t\right)\end{array}$
$+4\left[3\left(b-{x}_{EA}\right){h}^{\text{'}}\left(t\right)+\left(4{b}^{2}-6b{x}_{EA}+3{x}_{EA}^{2}\right){\alpha }^{\text{'}}\left(t\right)\right]\right\},$

herein, $b$ is the half-chord length of the airfoil, ${x}_{EA}=b{x}_{0}$ is streamwise position of the pitch axis measured from the leading edge, ${U}_{\infty }$, ${\rho }_{\infty }$ and ${M}_{\mathrm{\infty }}$ are the air speed, the air density, and the flight Mach number of the disturbed flow, respectively. $\lambda$ and $\gamma$ are correction aerodynamic factor and the isentropic gas coefficient, respectively.

Upon denoting ${x}_{1}=h/b$, ${x}_{2}=\alpha$, and the dimensionless time $\tau ={U}_{\infty }t/b$, Eq. (1) can be written as:

4
$\left\{\begin{array}{l}{\stackrel{˙}{x}}_{1}={x}_{3}\left(\tau \right),\\ {\stackrel{˙}{x}}_{2}={x}_{4}\left(\tau \right),\\ {\stackrel{˙}{x}}_{3}={a}_{1}^{\left(3\right)}{x}_{1}\left(\tau \right)+{a}_{2}^{\left(3\right)}{x}_{2}\left(\tau \right)+{a}_{3}^{\left(3\right)}{x}_{3}\left(\tau \right)+{a}_{4}^{\left(3\right)}{x}_{4}\left(\tau \right)+{\delta }_{A}{a}_{222}^{\left(3\right)}{x}_{2}^{3}\left(\tau \right)+{\delta }_{S}{a}_{222}^{\left(3\right)}{x}_{2}^{3}\left(\tau \right),\\ {\stackrel{˙}{x}}_{4}={a}_{1}^{\left(4\right)}{x}_{1}\left(\tau \right)+{a}_{2}^{\left(4\right)}{x}_{2}\left(\tau \right)+{a}_{3}^{\left(4\right)}{x}_{3}\left(\tau \right)+{a}_{4}^{\left(4\right)}{x}_{4}\left(\tau \right)+{\delta }_{A}{a}_{222}^{\left(4\right)}{x}_{2}^{3}\left(\tau \right)+{\delta }_{S}{a}_{222}^{\left(4\right)}{x}_{2}^{3}\left(\tau \right),\end{array}\right\$

where the coefficients can be found in Appendix, and the dots denote $d/d\tau$.

## 3. The stability and bifurcation analysis of system (4)

### 3.1. Stability of the initial equilibrium point

Obviously, $O\left(0,0,0,0\right)$ is an equilibrium point of Eq. (4). The Jacobian of the linearization system of system Eq. (4) evaluated at $O$ is:

5
$J\left(O\right)=\left(\begin{array}{l}0010\\ 0001\\ {a}_{1}^{\left(3\right)}{a}_{2}^{\left(3\right)}{a}_{3}^{\left(3\right)}{a}_{4}^{\left(3\right)}\\ {a}_{1}^{\left(4\right)}{a}_{2}^{\left(4\right)}{a}_{3}^{\left(4\right)}{a}_{4}^{\left(4\right)}\end{array}\right).$

The characteristic equation of the matrix (5) is:

6
${\omega }^{4}+p{\omega }^{3}+q{\omega }^{2}+r\omega +s=0,$

where:

7
$p=\frac{V\lambda \left[4+3{x}_{0}^{2}+6{x}_{0}\left({\chi }_{\alpha }-1\right)-6{\chi }_{\alpha }\right]+3{r}_{\alpha }^{2}\left(V\lambda +2M\mu {\zeta }_{\alpha }+2M\mu \varpi {\zeta }_{h}\right)}{3M\mu V\left({r}_{\alpha }^{2}-{\chi }_{\alpha }^{2}\right)},$
$q=\frac{3M\mu {r}_{\alpha }^{2}\left[M\mu \left(1+{\varpi }^{2}\right)+2{\zeta }_{\alpha }\left(V\lambda +2M\mu \varpi {\zeta }_{h}\right)\right]}{3{V}^{2}{M}^{2}{\mu }^{2}\left({r}_{\alpha }^{2}-{\chi }_{\alpha }^{2}\right)}$
$+\frac{V\lambda \left[V\left(\lambda +3M\mu -3M\mu {\chi }_{\alpha }\right)+8M\mu \varpi {\zeta }_{h}+6M\mu \varpi {x}_{0}^{2}{\zeta }_{h}-3M\mu {x}_{0}\left(V+4\varpi {\zeta }_{h}\right)\right]}{3{V}^{2}{M}^{2}{\mu }^{2}\left({r}_{\alpha }^{2}-{\chi }_{\alpha }^{2}\right)},$
$r=\frac{3{r}_{\alpha }^{2}\left(V\lambda +2M\mu {\varpi }^{2}{\zeta }_{\alpha }+2M\mu \varpi {\zeta }_{h}\right)+V\lambda \varpi \left[3\varpi {x}_{0}^{2}-6{x}_{0}\left(\varpi +V{\zeta }_{h}\right)+2\left(2\varpi +3V{\zeta }_{h}\right)\right]}{3M\mu {V}^{3}\left({r}_{\alpha }^{2}-{\chi }_{\alpha }^{2}\right)},$
$s=\frac{{\varpi }^{2}\left[M\mu {r}_{\alpha }^{2}+{V}^{2}\lambda \left(1-{x}_{0}\right)\right]}{M\mu {V}^{4}\left({r}_{\alpha }^{2}-{\chi }_{\alpha }^{2}\right)}.$

According to the Routh-Hruwitz criterion, the initial equilibrium solution $\left({x}_{1},{x}_{2},{x}_{3},{x}_{4}\right)=\left(0,0,0,0\right)$ is stable if the following conditions are satisfied:

8
$p>0,pq-r>0,s>0,r\left(pq-r\right)-{p}^{2}s>0.$

When conditions (8) are not satisfied, the initial equilibrium solution is unstable, and bifurcations may occur. First we give the following results on the eigenvalues of matrix (5).

Proposition 1. Zero is not an eigenvalue of matrix (5).

Proposition 2. There exist no two pairs of purely imaginary eigenvalues for matrix (5).

In fact, $r>0$ and $s>0$ due to all the parameters in the expressions of $r$, $s$ are positive and ${x}_{0}<1$. However suppose zero is an eigenvalue of matrix (5), then it follows, on the other hand, if there exist two pairs of purely imaginary eigenvalues for matrix (5), then it follows $r=0$. So Propositions 1 and 2 are obvious.

For the aeroelastic stability problems, the condition $sp/r-{p}^{2}/4>0$ should be satisfied. According to the above analysis, it is meaningful to consider the only one critical case: two eigenvalues with negative real parts and a pair of purely imaginary eigenvalues (this case corresponds to the limit cycle flutter). In this case, the following conditions are fulfilled:

9
$r\left(pq-r\right)-{p}^{2}s=0,sp/r-{p}^{2}/4>0.$

### 3.2. Normal form, universal unfolding and Hopf bifurcation

Choosing the following dimensionless values of parameters as [18]: $\mu =\text{100}$, ${\chi }_{\alpha }=\text{0.25}$, ${x}_{0}=\text{0.5}$, ${r}_{\alpha }=\text{0.5}$, $\varpi =\text{1.2}$, ${\zeta }_{\alpha }={\zeta }_{h}=\text{0}$ and $\lambda =\text{1}$, Eq. (4) can be rewritten as:

10
$\left(\begin{array}{l}{\stackrel{˙}{x}}_{1}\\ {\stackrel{˙}{x}}_{2}\\ {\stackrel{˙}{x}}_{3}\\ {\stackrel{˙}{x}}_{4}\end{array}\right)=A\left(\begin{array}{l}{x}_{1}\\ {x}_{2}\\ {x}_{3}\\ {x}_{4}\end{array}\right)+f{x}_{2}^{3},$

where:

$A=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ \frac{-1.92}{{V}^{2}}& \frac{50M-{V}^{2}}{150M{V}^{2}}& \frac{-1}{150M}& \frac{1}{900M}\\ \frac{1.92}{{V}^{2}}& \frac{-100M-{V}^{2}}{75M{V}^{2}}& \frac{-1}{75M}& \frac{-1.1}{45M}\end{array}\right),f=\left(\begin{array}{c}0\\ 0\\ \frac{-M\left(1+\gamma \right)}{1800}\\ \frac{-M\left(1+\gamma \right)}{900}\end{array}\right)+\left(\begin{array}{c}0\\ 0\\ \frac{B}{3{V}^{2}}\\ \frac{-B}{0.75{V}^{2}}\end{array}\right).$

Substituting the above dimensionless values of parameters and the Eq. (7) into the first formula of (9), the speed on the flutter instability boundary is:

11
${V}_{F}=\sqrt{\frac{7.850112004×1{0}^{11}{M}^{2}}{-3.90656×1{0}^{8}+1.5859199×1{0}^{10}M}}.$

The relation of the flutter speed ${V}_{F}$ and flight Mach number $M$ is depicted in Fig. 2. The figure reveals that the increase of flight Mach number can result in the increase of the flutter safety.

Fig. 2The flutter speed VF vs. flight Mach number M

In the following, we assume $O$ is a bifurcation point, i.e. the conditions (9) hold. Substituting $V={V}_{F}$ into Eq. (10), we get:

12
$\left(\begin{array}{l}{\stackrel{˙}{x}}_{1}\\ {\stackrel{˙}{x}}_{2}\\ {\stackrel{˙}{x}}_{3}\\ {\stackrel{˙}{x}}_{4}\end{array}\right)=\stackrel{-}{A}\left(\begin{array}{l}{x}_{1}\\ {x}_{2}\\ {x}_{3}\\ {x}_{4}\end{array}\right)+\stackrel{-}{f}{x}_{2}^{3},$

where:

$\stackrel{-}{A}=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ \frac{-1.92}{{V}_{F}^{2}}& \frac{50M-{V}_{F}^{2}}{150M{V}_{F}^{2}}& \frac{-1}{150M}& \frac{1}{900M}\\ \frac{1.92}{{V}_{F}^{2}}& \frac{-100M-{V}_{F}^{2}}{75M{V}_{F}^{2}}& \frac{-1}{75M}& \frac{-1.1}{45M}\end{array}\right),\stackrel{-}{f}=\left(\begin{array}{c}0\\ 0\\ \frac{-M\left(1+\gamma \right)}{1800}\\ \frac{-M\left(1+\gamma \right)}{900}\end{array}\right)+\left(\begin{array}{c}0\\ 0\\ \frac{B}{3{{V}_{F}}^{2}}\\ \frac{-B}{0.75{{V}_{F}}^{2}}\end{array}\right).$

For Eq. (12), there is a coordinate transformation $x=Py\text{,}$$x=\left({x}_{1},{x}_{2},{x}_{3},{x}_{4}{\right)}^{T}\text{,}$$y=\left({y}_{1},{y}_{2},{y}_{3},{y}_{4}{\right)}^{T}$ so that:

13
$\stackrel{˙}{y}=Jy+g{\left(\sum _{j=1}^{4}{P}_{2j}{y}_{j}\right)}^{3},$

where $J={P}^{-1}\stackrel{-}{A}P=\left(\begin{array}{l}0-\omega 00\\ \omega 000\\ 00-\alpha -\beta \\ 00\beta -\alpha \end{array}\right)$, $±\mathrm{i}\omega$ and $-\alpha ±\mathrm{i}\beta$ are the eigenvalues of the matrix $\stackrel{-}{A}$, and $g=\left({g}_{1},{g}_{2},{g}_{3},{g}_{4}{\right)}^{T}={P}^{-1}\stackrel{-}{f}$.

Using the Maple program in reference [19], the normal form of Eq. (13) is calculated as:

14
$\left(\begin{array}{l}{\stackrel{˙}{z}}_{1}\\ {\stackrel{˙}{z}}_{2}\end{array}\right)=\left(\begin{array}{ll}0& -\omega \\ \omega & 0\end{array}\right)\left(\begin{array}{l}{z}_{1}\\ {z}_{2}\end{array}\right)+\left({z}_{1}^{2}+{z}_{2}^{2}\right)\left\{{\delta }_{1}\left(\begin{array}{l}{z}_{1}\\ {z}_{2}\end{array}\right)+{\delta }_{2}\left(\begin{array}{l}-{z}_{2}\\ {z}_{1}\end{array}\right)}\right\,$

where:

15
${\delta }_{1}=\frac{3}{8}\left({g}_{1}{P}_{21}{P}_{22}^{2}+{g}_{1}{P}_{21}^{3}+{g}_{2}{P}_{22}{P}_{21}^{2}+{g}_{2}{P}_{22}^{3}\right),$
${\delta }_{2}=\frac{-3}{8}\left({g}_{1}{P}_{22}{P}_{21}^{2}+{g}_{1}{P}_{22}^{3}-{g}_{2}{P}_{21}{P}_{22}^{2}-{g}_{2}{P}_{21}^{3}\right).$

The coordinate transformation $y=Tz$ is omitted here (for details, see the Maple program in reference [19]).

If ${\delta }_{1}\ne 0$, the universal unfolding of Eq. (14) can be written as [20]:

16
$\left(\begin{array}{l}{\stackrel{˙}{z}}_{1}\\ {\stackrel{˙}{z}}_{2}\end{array}\right)=\left(\begin{array}{ll}\nu & -\omega \\ \omega & \nu \end{array}\right)\left(\begin{array}{l}{z}_{1}\\ {z}_{2}\end{array}\right)+\left({z}_{1}^{2}+{z}_{2}^{2}\right)\left\{{\delta }_{1}\left(\begin{array}{l}{z}_{1}\\ {z}_{2}\end{array}\right)+{\delta }_{2}\left(\begin{array}{l}-{z}_{2}\\ {z}_{1}\end{array}\right)}\right\.$

Because Eq. (16) is the universal unfolding of Eq. (14), Eq. (16) bears all the dynamic behavior of any small perturbed system of Eq. (14) in the vicinity of $z={\left({z}_{1},{z}_{2}\right)}^{T}=0$. Furthermore, using the transformations $P$ and $T$, we can obtain all the dynamic behavior of any small perturbed system of Eq. (12) in the vicinity of $O$.

Let ${z}_{1}=r\mathrm{c}\mathrm{o}\mathrm{s}\theta$ and ${z}_{2}=r\mathrm{s}\mathrm{i}\mathrm{n}\theta$, then the polar coordinate form of Eq. (16) is:

17
$\left\{\begin{array}{l}\stackrel{˙}{r}=r\left(\nu +{\delta }_{1}{r}^{2}\right),\\ \stackrel{˙}{\theta }=\omega +{\delta }_{2}{r}^{2}.\end{array}\right\$

Now we derive the expression of the unfolding parameter $\nu$ in Eq. (17).

Consider one of the perturbed systems of Eq. (12) as follows. Without loss of generality, using the parameter transformations $V={V}_{F}+\eta$ and the state variable transformation $x=Py$, one may transform Eq. (12) into a new system as follows:

18
$\stackrel{˙}{y}={P}^{-1}{\stackrel{-}{A}}_{\eta }Py+{P}^{-1}{\stackrel{-}{f}}_{\eta }\left(\sum _{j=1}^{4}{P}_{2j}{y}_{j}{\right)}^{3},$

where:

${\stackrel{-}{A}}_{\eta }=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ \frac{-1.92}{{\left({V}_{F}+\eta \right)}^{2}}& \frac{50M-{\left({V}_{F}+\eta \right)}^{2}}{150M{\left({V}_{F}+\eta \right)}^{2}}& \frac{-1}{150M}& \frac{1}{900M}\\ \frac{1.92}{{\left({V}_{F}+\eta \right)}^{2}}& \frac{-100M-{\left({V}_{F}+\eta \right)}^{2}}{75M{\left({V}_{F}+\eta \right)}^{2}}& \frac{-1}{75M}& \frac{-1.1}{45M}\end{array}\right),$
${\stackrel{-}{f}}_{\eta }=\left(\begin{array}{c}0\\ 0\\ \frac{-M\left(1+\gamma \right)}{1800}\\ \frac{-M\left(1+\gamma \right)}{900}\end{array}\right)+\left(\begin{array}{c}0\\ 0\\ \frac{B}{3{\left({V}_{F}+\eta \right)}^{2}}\\ \frac{-B}{0.75{\left({V}_{F}+\eta \right)}^{2}}\end{array}\right).$

Let $L=\left({L}_{1},{L}_{2},{L}_{3},{L}_{4}{\right)}^{T}={P}^{-1}{\stackrel{-}{A}}_{\eta }Py$. According to reference [21], in Eq. (17) the unfolding parameter $\nu =\frac{1}{2}\eta ×{\left(\frac{{\partial }^{2}{L}_{1}}{\partial {y}_{1}\partial \eta }+\frac{{\partial }^{2}{L}_{2}}{\partial {y}_{2}\partial \eta }\right)}_{{y}_{1}={y}_{2}=\eta =0}$. We can obtain the conclusions for Eqs. (16)-(18) as follows:

Case 1: ${\delta }_{1}<\text{0}$.

If $\nu <\text{0}$, (0,0) is a stable focus of Eq. (16), so (0,0,0,0) is also a stable focus of Eq. (18); if $\nu =\text{0}$, because the solutions of Eq. (17), ${r}^{2}=-1/2\left({\delta }_{1}t+c\right)$, are attracted to zero when $t\to +\mathrm{\infty }$, (0,0,0,0) is an asymptotically stable, but non-hyperbolic equilibrium point of Eq. (18); if $\nu >\text{0}$, for Eq. (16), (0,0) is an unstable focus and there is a stable LCO , ${z}_{1}^{2}+{z}_{2}^{2}=-v/{\delta }_{1}$, so for Eq. (18) (0,0,0,0) is an unstable focus and there is a stable LCO as well. Supercritical Hopf bifurcation takes place in this case. In practice, because the amplitude of the LCO builds up from zero gradually as $\nu$ increases, it doesn’t bring immediately failure of the structure. So the flutter instability in this case is benign.

Case 2: ${\delta }_{1}>\text{0}$.

For Eq. (18), if $\nu <\text{0}$, (0,0,0,0) is a local stable focus, and there is an unstable LCO; if $\nu =\text{0}$, (0,0,0,0) is an unstable non-hyperbolic equilibrium point; if $\nu >0$, (0,0,0,0) is an unstable focus. Subcritical Hopf bifurcation occurs in this case. In practice, when $\nu >\text{0}$, the response with any small initial conditions will move away from the trivial point immediately and the flutter is violent. So the flutter instability in this case is catastrophic.

## 4. The effect of the structural and aerodynamic parameters on character of the flutter instability

One concludes that the sign of ${\delta }_{1}$ decides the type of Hopf bifurcation or the character of flutter instability, i.e. benign or catastrophic. ${\delta }_{1}$ vs. $M$ as a function of the normalized nonlinear stiffness coefficient $B$ when $\gamma =\text{1.4}$ are shown in Fig. 3 (curves are drawn pointwisely by using the first formula of Eq. (15)).

Fig. 3The relations of δ1 and flight Mach number M as a function of the structural nonlinearity B

The Fig. 3 indicates that the aerodynamic nonlinearity, representation $M$, and the structural nonlinearity, representation $B$, can affect the sign of ${\delta }_{1}$.

For $B=-\text{50}$, $-\text{30}$ and $\text{0}$, the whole curves ${\delta }_{1}-M$ lie above the $M$ axis. So the flutter takes place always as a subcritical Hopf bifurcation for any value of $M$, and the flutter instability is catastrophic. For $B=\text{30}$ and 50, the curves ${\delta }_{1}-M$ intersect the $M$ axis at the critical values. We define the critical value as ${M}_{c}$, at where the sign of ${\delta }_{1}$ changes from negative to positive. This implies that with the increase of the flight Mach number $M$, the flutter instability changes from benign to catastrophic at last. In other respect, the increase of hard structural nonlinearity ($B>0$) results in the increase of ${M}_{c}$ and the flight safety, and for hard structural nonlinearity the decrease of flight Mach number $M$ can increase the flight safety as well.

## 5. Examples, numerical simulations, and the effect of the parameters on character of the flutter instability

In the following, two cases with the flight Mach number $M=\text{4}$ and $M=\text{10}$ are considered and the flutter speeds are solved from Eq. (11) as ${V}_{1F}=\text{14.11460254}$ for $M=\text{4,}$ and ${V}_{2F}=\text{22.27577602}$ for $M=\text{10}$, respectively.

Take $M=\text{4}$ for example. In this case, the parameters in Eq. (16) were determined as:

${\delta }_{1}=4.519094258+4.519094258\gamma -3.278155804B,$
${\delta }_{2}=5.495174328+5.495174328\gamma -3.181798633B,$
$\nu =0.006589997791\eta .$

Fig. 4The relation of δ1, B and γ for M=4

Fig. 5Supercritical Hopf bifurcation for γ=1.4, B=100, M=4: a) trajectory starting from x0T=(0.5,0.28,0,0) converges to O for V=14, b) the transient time history of h/b for V=14.3 and x0T=(0.0001,0.0001,0,0), c) the transient time history of α for V=14.3 and x0T=(0.0001,0.0001,0,0), d) trajectory starting from x0T=(0.5,0.28,0,0) converges to a LCO for V=14.3

a)

b)

c)

d)

The relation of ${\delta }_{1}$, $B$ and $\gamma$ is depicted in Fig. 4. The character of flutter instability is benign above the slope line, but it is catastrophic below the slope line. The decrease of isentropic gas coefficient $\gamma$ results in the increase of the flight safety, and the increase of hard structural nonlinearity ($B>\text{0}$) also does.

Fig. 6Subcritical Hopf bifurcation for γ=1.4, B=2.5, M=4: a) trajectory starting from x0T=(0.04,0.0001,0,0) converges to a LCO immediately for V=14, b) trajectory starting from x0T=(0.0002,0.00015,0,0) converges to a LCO for V=14.3

a)

b)

Fig. 7Supercritical Hopf bifurcation for γ=1.4, B=21, M=10: a) trajectory starting from x0T=(0.8,0.5,0,0) converges to O for V=22.1, b) the transient time history of h/b for V=22.5 and x0T=(0.001,0.005,0,0), c) the transient time history of α for V=22.5 and x0T=(0.8,0.7,0,0)

a)

b)

c)

To verify these theoretical results, numerical simulations using the Runge-Kutta algorithm to the Eq. (10) for $M=\text{4}$ were carried out. Notice that ${x}_{1}=h/b$, ${x}_{2}=\alpha$. When $B=\text{100}$ and $\gamma =\text{1.4}$, ${\delta }_{1}<\text{0}$. Supercritical Hopf bifurcation occurs, as shown in Fig. 5. Because the amplitude of the LCO is small no matter the initial distance is small or large, the flutter instability is benign. When $B=\text{2.5}$ and $\gamma =\text{1.4}$, ${\delta }_{1}>\text{0}$. Subcritical Hopf bifurcation occurs, as shown in Fig. 6. The amplitudes of the two LCOs are large, so the flutter instability yields catastrophic failure of the structure.

For the case with the flight Mach number $M=\text{10}$, the flutter speed ${V}_{2F}=\text{22.27577602}$, ${\delta }_{1}=\text{27.18679908}+\text{27.18679908}\gamma -\text{3.159223268}B\text{,}\text{}\text{}$${\delta }_{2}=\text{51.27974352}+\text{51.27974352}\gamma -\text{5.056781381}B$, and $\nu =\text{0.00227926149}\eta$. The relation of ${\delta }_{1}$, $B$ and $\gamma$ is similar to what is depicted in Fig. 4. Numerical simulations for the Eq. (10) when $M=\text{10}$ were carried out. When $B=\text{21}$ and $\gamma =\text{1.4}$, ${\delta }_{1}<\text{0}$, Fig. 7 indicates that the bifurcation is supercritical and the flutter instability is benign. Whilst when $B=-\text{10}$ and $\gamma =\text{1.4,}$${\delta }_{1}>\text{0}$ Fig. 8 indicates that the subcritical Hopf bifurcation takes place and the flutter instability is catastrophic.

Fig. 8Subcritical Hopf bifurcation for γ=1.4, B=-10, M= 10: a) trajectory starting from x0T=(0.001,0,0,0) converges to O for V=22.1, b) trajectory starting from x0T=(0.00002,0.00001,0,0) moves away from O immediately for V=22.5

a)

b)

## 6. Conclusions

This paper deals with the flutter of a two-dimensional lifting surfaces in a supersonic flow field. This Eq. (4) undergoes only Hopf bifurcations in the neighborhood of the trivial equilibrium point if this system bifurcates at this point. The normal form and universal unfolding are calculated by using the maple program and the bifurcation theory. The supercritical Hopf bifurcation is benign, because it brings less flutter to the airfoil. Comparably the subcritical Hopf bifurcation is catastrophic, because it renders a harmful violent flutter of the airfoil.

In addition, we study the influence of the structural and aerodynamic parameters, representations $B\text{,}$$M$ and $\gamma \text{,}$ on the character of flutter instability. The influence of these parameters can be described as follows:

1) The soft structural nonlinearities ($B<\text{0}$) yield failure of the structure (i.e. catastrophic flutter), and yet the hard structural nonlinearities ($B>\text{0}$) perhaps not.

2) The increase of hard structural nonlinearity results in the increase of the airfoil safety, and for hard structural nonlinearity the decrease of the flight Mach number $M$ and isentropic gas coefficient $\gamma$ also do.

#### References

• Zhang Q. C., Liu H. Y., Ren A. D. The study of limit cycle flutter for airfoil with nonlinearity. Acta Aerodynamica Sinica, Vol. 22, Issue 3, 2004, p. 332-336, (in Chinese).
• Yang X. L., Sethna P. R. Local and global bifurcations in parametrically excited vibrations of nearly square plates. International Journal of Non-Linear Mechanics, Vol. 26, Issue 2, 1991, p. 199-220.
• Zhang Q. C., Liu H. Y., Ren A. D. Local bifurcation for airfoil with cubic nonlinearities. Journal of Tianjin University, Vol. 37, Issue 11, 2004, p. 970-974, (in Chinese).
• Breitbach E. Effect of structural nonlinearities on aircraft vibration and flutter. The 45th Structures and Materials AGARD Panel Meeting, Agard report 665, 1977.
• Librescu L. Elastostatics and Kinetics of Anisotropic and Heterogeneous Shell-Type Structures, Aeroelastic Stability of Anisotropic Multilayered Thin Panels. Mechanics of Elastic Stability. Noordhoff International Publishing, Leyden, 1975, p. 53-63, p. 106-158, 543-550.
• Dowell E. H., Ilgamov M. Studies in Non-Linear Aeroelasticity. Springer-Verlag, New York, 1988.
• Raghothama A., Narayanan S. Non-linear dynamics of a two-dimensional airfoil by incremental harmonic balance method. Journal of Sound and Vibration, Vol. 226, Issue 3, 1999, p. 493-517.
• Shahrzad P., Mahzoon M. Limit cycle flutter of airfoils in steady and unsteady flows. Journal of Sound and Vibration, Vol. 256, Issue 2, 2002, p. 213-225.
• Liu L., Wong Y. S., Lee B. H. K. Non-linear aeroelastic analysis using the point transformation method, Part 1: freeplay model. Journal of Sound and Vibration, Vol. 253, Issue 2, 2002, p. 447-469.
• Zhao D. M., Zhang Q. C. Bifurcation and chaos analysis for aeroelastic airfoil with freeplay structural nonlinearity in pitch. Chinese Physics B, Vol. 19, Issue 3, 2010, p. 1-10.
• Yang Y. R. KBM method of analyzing limit cycle flutter of a wing with an external store and comparison with a wind-tunnel test. Journal of Sound and Vibration, Vol. 187, Issue 2, 1995, p. 271-280.
• Kim S. H., Lee I. Aeroelastic analysis of a flexible airfoil with a freeplay non-linearity. Journal of Sound and Vibration, Vol. 193, Issue 4, 1996, p. 823-846.
• Price S. J., Alighanbari H., Lee B. H. K. The aeroelastic response of a two-dimensional airfoil with bilinear and cubic structural nonlinearities. Journal of Fluids and Structures, Vol. 9, Issue 2, 1995, p. 175-193.
• Abbas L. K., Chen Q., Donnell K., Valentine D., Marzocca P. Numerical studies of a non-linear aeroelastic system with plunging and pitching freeplays in supersonic/hypersonic regimes. Aerospace Science and Technology, Vol. 11, Issue 5, 2007, p. 405-418.
• Liu L., Wong Y. S., Lee B. H. K. Application of the center manifold theory in non-linear aeroelasticity. Journal of Sound and Vibration, Vol. 234, Issue 4, 2000, p. 641-659.
• Ding Q., Wang D. L. The flutter of an airfoil with cubic structural and aerodynamic non-linearities. Aerospace Science and Technology, Vol. 10, Issue 5, 2006, p. 427-434.
• Lee B. H. K., Liu L., Chung K. W. Airfoil motion in subsonic flow with strong cubic nonlinear restoring forces. Journal of Sound and Vibration, Vol. 281, Issue 3, 2005, p. 699-717.
• Librescu L., Chiocchia G., Marzocca P. Implications of cubic physical/aerodynamic non-linearities on the character of the flutter instability boundary. International Journal of Non-Linear Mechanics, Vol. 38, Issue 2, 2003, p. 173-199.
• Bi Q. S., Yu P. Symbolic computation of normal forms for semi-simple cases. Journal of Computational and Applied Mathematics, Vol. 102, Issue 2, 1999, p. 195-220.
• Lu Q. S. Advanced Series in Nonlinear Science Bifurcation and Singularity. Shanghai Scientific and Technological Education Publishing House, Shanghai, 1995, (in Chinese).
• Yu P. Computation of normal forms via a perturbation technique. Journal of Sound and Vibration, Vol. 211, Issue 1, 1998, p. 19-38.