Analysing Non-Linear Flutter Vibrations Using System Dynamic Approach

Abstract

The objective of this work is to investigate the dynamic behaviour of aero-elastic vibrations in the presence of non-linear stiffness such as free-play mechanism, softening or hardening stiffness. A closed loop dynamic system is proposed to represent the phenomenon of flow-structure interaction. In this approach a transfer function for generating the aerodynamic forces based on the structural response is constructed in the feedback loop of the dynamic system with the aid of Padé rational function. The effects of the non-linear factors therefore can be included conveniently in time domain simulation and the stability of limit cycle oscillations (LCO) can be analysed accurately.

1. Introduction

Aero-elasticity is a multi-physics discipline that involves the loads of aerodynamics, elastic and inertial generated by the motion of structure. One of the most important phenomena in this field is flutter regarding to its harmful effect to the structure. This flow-induced vibration under certain conditions can be self-excited and divergently unstable. In aerospace industry the boundary of flutter instability is usually determined by V-g method, a computational technique in frequency domain based on the balance of energy of the oscillating wing and the flow by maintaining a harmonic function of the aero-elastic response. In this method explained remarkably by Stanciu et al. [1], the critical velocity of the flow is determined by solving the complex eigen value problem of the aero-elastic system. Although the method can give accurate results it is only effective for linear cases. However in real aircraft there are non-linear factors of structure such as free-play, hysteretic and large deformation that need to be taken into account. Trickey et al. [2] observed that certain cases in regards to excessive LCO found in some Boeing and Airbus aircrafts and stated that the characterization and explanation of this non-linear vibration were important for fatigue and maintenance issues. In their research, the methods of non-linear dynamics were developed for these purposes and they proposed a novel system-identification technique to generate an approximation of LCO to be used for online monitoring of dynamic behaviour close to bifurcation condition [8].

Pereira et al. [3] showed an example of LCO due to the existence of hardening nonlinearity of wing stiffness in pitching of F-16 aircraft that caused persistent aero-elastic problems. Therefore the knowledge and comprehension of non-linear aero-elasticity are of increasing importance in aircraft design. In their work an investigation on the combined influence of hardening and free-play nonlinearities on the bifurcation response was carried out.

A document regarding the missing of MH370 Boeing 777 is also concerned about the phenomena of aero-elasticity. The failure analysis of the right-side flaperon that was found in French territory’s Reunion Island [4] on 2015, reported that flutter (LCO) caused to repetitive loading which in turn imparted stress fatigue in the primary aluminum alloy attachment components.

This work has an intention to simulate numerically the interaction of flow-structure as a dynamic system by arranging the structure part as a principle plant and the aerodynamic part as a feedback loop subsystem. The analysis of structural response in time domain enables to insert the nonlinearities conveniently. The part of structures is reconstructed in a form of block-diagram representing a dynamic system where the inertial loads are expressed explicitly as a result of the elastic loads and frictions generated in the progressing structure response due to external aerodynamic excitations, while the part of aerodynamics is arranged as a feedback-loop transfer function activated by the structural response. For this purpose, the unsteady aerodynamic forces calculated by using singularity method in frequency domain have to be converted to Laplace variable s by using Padé’s approximation rational function. Botez et al. [5] in conducting flutter analysis of CL-604 Bombardier, used a least-squares technique utilizing certain number of lagging-terms, and the approximation showed the best aerodynamic forces conversion from frequency into Laplace domain in terms of execution time and precision.

In analysing LCO on mechanical system in general where there is involvement of various physical parameters such as non-linear stiffness, hysteretic and free-play, and more specifically the influence of damping to the stability of the oscillations, Sinou and Jézéquel [6] proposed to employ a two-degree-of-freedom model for the sake of simplicity. With the same spirit, in this study we use a pitch-plunge two dimensional wing-section model in analysing the effects of structural nonlinearities on a binary classical flutter.

2. Stability analysis of aero-elastic system

2.1 Description of the two-degree-of-freedom model

Figure 1 shows schematically a two-degree-of-freedom pitch-plunge fluttering aerofoil model. The support system consisting of axial and rotational springs are attached to a rigid aerofoil on a point so-called elastic axis. These two flexible supports restrict the motions of the aerofoil with the exception in the two modes of translation and rotation. For the case of zero damping, the equations of motion of the aerofoil subjected a uniform flow can be written as:

$$\left[\begin{array}{cc}m& {\mathit{mx}}_{\alpha }b\\ {\mathit{mx}}_{\alpha }b& I_g+m{\left(x_{\alpha }b\right)}^2\end{array}\right]\left\{\begin{array}{c}\ddot{h}\\ \ddot{\alpha }\end{array}\right\}+\left[\begin{array}{cc}{k}_h& 0\\ 0& k_{\alpha }\end{array}\right]\left\{\begin{array}{c}h\\ \alpha \end{array}\right\}=\left\{\begin{array}{c}-L\\ \mathit{Lec}+M_{\mathit{ac}}\end{array}\right\}=\left\{F_{\mathit{\text{aero}}}\right\}$$ E1

where h and α are the degree of freedom in plunging and pitching, respectively, as explained in Figure 1, while k_h and k_α are the spring stiffness in translation and rotation, respectively.

Figure 1.

Two-degree-of-freedom aero-elastic model.

2.2 Unsteady aerodynamic model

Theodorsen’s unsteady aerodynamic model as explained by Brunton and Rowly [7] excellently is used in this work. This method analyses the motions of the aerofoil in frequency domain and it assumes that the amplitudes are small. In the analysis the aerofoil is considered thin, the flow is inviscid incompressible with no separation or intrusion.

In this method the frequencies of the harmonic oscillating motions are considered relatively slow therefore the transversal and rotational velocities of the aerofoil contribute as an additional angle of attack to the total lift. As a result the quasi-steady of the lift coefficient $C_L^{\mathit{QS}}$ can be expressed proportional to the total angle of attack:

$$C_L^{\mathit{QS}}=\frac{\partial C_l}{\partial \alpha }\left(\alpha +\frac{\dot{h}}{U_{\infty }}+b\left(\frac{1}{2}-a\right)\frac{\dot{\alpha }}{U_{\infty }}\right)$$ E2

In thin aerofoil theory the lift gradient can be considered equals to 2π, the vortex singularity is located at the aerodynamic center (a quarter of the chord from the leading-edge) and the downwash velocity is focused at three quarter of the chord.

The aerodynamic loading consisted of lift and pitching moment can be presented as:

$$L=\frac{1}{2}\rho U^{2}clC\left(k\right)C_L^{\mathit{QS}}$$ E3

$$M_{\mathit{ac}}=\frac{1}{2}\rho U^2{c}^{2}lC\left(k\right)C_L^{\mathit{QS}}$$ E4

where c and l are the chord and the span of the 2D wing, respectively. C(k) is Theodorsen’s function showing that there is a phase difference between aerodynamic loading and wing section’s motion. The parameter k called reduced (nondimensional) frequency is defined as ωb/U_∞.

2.3 Flutter stability boundary

V-g method based on the balance of energy between the flow and the motion of the structure is used to determine the flutter boundary for a linear aero-elastic system. The analysis is conducted in frequency domain where harmonic motions in pitching and plunging are imposed to the dynamic response of the structure to represent the state of the aero-elastic system in the stability boundary. In order to maintain harmonic motions of the structure, a virtual structural damping g is inserted to the system hence the equation of motions adopting from Eq. (1) becomes:

$$\left[\begin{array}{cc}m& {\mathit{mx}}_{\alpha }b\\ {\mathit{mx}}_{\alpha }b& I_g+m{\left(x_{\alpha }b\right)}^2\end{array}\right]\left\{\begin{array}{c}\ddot{h}\\ \ddot{\alpha }\end{array}\right\}+\left(1+\mathit{ig}\right)\left[\begin{array}{cc}{k}_h& 0\\ 0& k_{\alpha }\end{array}\right]\left\{\begin{array}{c}h\\ \alpha \end{array}\right\}=\left\{F_{\mathit{\text{aero}}}\right\}$$ E5

By imposing a harmonic functions to the motions in both plunging and pitching as shown in Eq. (6):

$$\left\{\begin{array}{c}h\\ \alpha \end{array}\right\}=\left\{\begin{array}{c}\widehat{h}\\ \widehat{\alpha }\end{array}\right\}{e}^{\mathit{i\omega t}}$$ E6

Henceforth, the aerodynamic loading can be represented as the expression in Eq. (7):

$$\left\{F_{\mathit{\text{aero}}}\right\}=\left\{\begin{array}{c}-L\\ \mathit{Lec}+M_{\mathit{ac}}\end{array}\right\}=\frac{1}{2}\rho U^2\left[Q\left(\mathit{ik}\right)\right]\left\{\begin{array}{c}h\\ \alpha \end{array}\right\}$$ E7

where the generalized aerodynamic matrix Q(ik) in complex form can be written as:

$$\left[Q\left(\mathit{ik}\right)\right]=\mathit{cl}C\left(k\right)\times \left(\left[\begin{array}{cc}0& \left(-\frac{\partial C_l}{\partial \alpha }\right)\\ 0& \left(2\frac{\partial C_l}{\partial \alpha }\mathit{be}+2\frac{\partial C_m}{\partial \alpha }b\right)\end{array}\right]+i\left[\begin{array}{cc}\left(-\frac{\partial C_l}{\partial \alpha }\frac{k}{b}\right)& 0\\ \left(2\mathit{ek}\frac{\partial C_l}{\partial \alpha }+2k\frac{\partial C_m}{\partial \alpha }\right)& 0\end{array}\right]\right)$$ E8

subsequently the equations of motion showed in Eq. (5) will lead to the solution of the eigen values problem as presented below:

$${\left[\begin{array}{cc}{k}_h& 0\\ 0& k_{\alpha }\end{array}\right]}^{-1}\left(\left[\begin{array}{cc}m& {\mathit{mx}}_{\alpha }b\\ {\mathit{mx}}_{\alpha }b& I_g+m{\left(x_{\alpha }b\right)}^2\end{array}\right]+\frac{\rho b^2}{2k^2}\left[Q\left(\mathit{ik}\right)\right]\right)\left\{\begin{array}{c}\widehat{h}\\ \widehat{\alpha }\end{array}\right\}=\frac{\left(1+\mathit{ig}\right)}{{\omega }^2}\left\{\begin{array}{c}\widehat{h}\\ \widehat{\alpha }\end{array}\right\}$$ E9

Elaborating a certain range of reduced frequency k into Eq. (9) will give as a result, a range of complex eigenvalues for both modes of motion. The real parts relate to the natural frequencies while the imaginary parts to the artificial structural damping of the aero-elastic system for certain values of flow velocity, associated with the reduced frequency k (Table 1).

Notation	Description	Value
a	Relative distance to half chord from midchord	−0.5
b	Half chord	0.05 m
c	Chord	0.1 m
l	Span of the wing model	0.28 m
e	Relative distance from E.A. to a.c.	0.0
xα	Relative distance from c.g. to mid-chord in b	0.34
m	Plunging mass	1.35 kg
mp	Pitching mass	0.93 kg
Ip	Mass moment of inertia in pitching	7.741e^−4 kg m^2
kh	Axial stiffness in plunging	906.81 N/m
kα	Rotational stiffness in pitching	11,445 Nm
E.A.	Elastic axis
a.c.	Aerodynamic center
c.g.	Centre of gravity

Table 1.

Value of physical parameters.

An aero-elastic flutter model with NACA0015 wing section is fabricated and installed in the wind-tunnel in vertical position as shown in Figure 2a and b. The each end of the wing model is mounted on a support system consisted of a pair of steel cantilever beams to allow the side-slipping translation motion and a warm spring for the yawing rotational motion.

Figure 2.

Aero-elastic wind-tunnel model, NACA 0015 aerofoil, two-degree-of-freedom system in translation and rotation, in the cross-section of 30 cm × 30 cm test-section.

The solution of the eigenvalue problem in Eq. (9) yields two curves of natural frequencies and two artificial structural damping depending to the flow speed as presented in Figures 3 and 4, respectively. The critical speed is defined where one of the artificial structural damping curves intercepts the real actual structural damping of the structure.

Figure 3.

Natural frequencies in translation and rotation modes (in Hertz) versus flow velocity (in m/s).

Figure 4.

Artificial structural damping parameters g in translation and rotation modes versus flow velocity (in m/s). The solid horizontal red colour line is the value of structural damping of steel (the material of the both springs).

Figure 4 shows the flutter boundary of this linear aero-elastic model where the critical flow speed is around 15 m/s. The negative values of the damping curves indicate that an amount of energy has to be supplied to the system to have a harmonic response and for the positive ones the energy has to be dissipated such to maintain the harmonic stable response. In Figure 3, the frequencies of two aero-elastic modes, i.e. translation and rotation approach each other as the flow speed is getting close to the critical speed of flutter boundary. The phenomenon so-called internal resonance shows that there is an interchange of energy between the two modes of vibration.

3. Non-linear behaviour of the aero-elastic system

3.1 System dynamic approach

Consider a structure with M, C, K as the matrices of mass, damping and stiffness respectively, subjected to an external loads F(t). The dynamic response x(t) basically can be presented in an arrangement of block diagrams based on the equation of motion by showing explicitly the inertial internal loads:

$$M\ddot{x}=F\left(t\right)-C\dot{x}-Kx$$ E10

The above expression can be considered as a junction with the input of F(t) which consists of a disturbance and aerodynamic forces encountered by closed loop feedback signals of C and K, brings forth the output of inertial load signals. Figure 5 explains the block diagram of this mechanical system.

Figure 5.

Block diagram for simulation of a structure subjected to an external load and some initial conditions.

To be able to model the effects of structural non-linearity accurately and to simulate conveniently in time domain, the aero-elastic system needs to be represented as a system dynamic model. In this system the aerodynamic forces are generated by the lifting surface as a result of the structural temporal response in a closed loop form. For this purpose the unsteady-aerodynamic forces calculated in frequency domain based on harmonic motions of the natural modes as expressed in Eq. (8), are now converted in Laplace variable using Padé rational function approximation as shown in Eq. (11).

The lagging term parameters β_j, for j = 1, n, are real and chosen less than 1, where the values and the numbers are determined to optimize the approximated curves. The matrices A₀, A₁ and A₂ in real values are estimated with curve-fitting using least square technique in complex plan to approximate the values of the aerodynamic forces calculated in frequency domain.

$$\left[Q\left(\mathit{ik}\right)\right]\approx \left[A_0\right]+\left[A_1\right]\left(\frac{\mathit{bs}}{U_{\infty }}\right)+\left[A_2\right]{\left(\frac{\mathit{bs}}{U_{\infty }}\right)}^2+\sum _{j=1}^n\frac{\left[A_{j+2}\right]s}{s+{\beta }_j\frac{U_{\infty }}{b}}$$ E11

Figure 6A and B explains the curve fitting of the generalised aerodynamic forces in matrix Q11 and Q21 calculated in a range of frequency in discretized data with the approximation rational function of Padé showed in solid lines.

Figure 6.

(A) Matrix Q11 signifying the generalized aerodynamic forces where the points are values of calculated in frequency domain for a certain range of reduced frequency k, and the solid line is the approximated curve in Laplace variable s by using four parameters of lagging term β_j of 01, 0.2, 0.3 and 0.4. (B) Matrix Q21 signifying the generalized aerodynamic forces where the points are values of calculated in frequency domain for a certain range of reduced frequency k, and the solid line is the approximated curve in Laplace variable s by using four parameters of lagging term β_j of 01, 0.2, 0.3 and 0.4.

The other matrix of the generalised aerodynamic forces Q12 and Q22 are zero as the consequence of the location of the elastic axis coincides with the aerodynamic center of the wing section.

Based on the dynamic response of structure in terms of x(t), d/dt x(t) and d2/dt2 x(t) the aerodynamic forces in the forms of lift and pitching moment can be calculated by using a transfer function constructed using Eq. (10). Figure 7 shows the arrangement of the block diagrams to express the aerodynamic transfer function.

Figure 7.

Aerodynamic transfer function.

The interaction of flow and structure can be simulated by arranging the block diagrams of structure showed in Figure 5 as a subsystem, coupled with the aerodynamic forces showed in Figure 7 where the calculation is conducted in time domain for each step of time discrete, simultaneously. The block diagrams of the plant representing the structure and the aerodynamic forces subsystem as a feedback loop have to be arranged such that all the processes are enhanced in integral operations to ensure the minimum numerical errors and the convergence of the solutions.

By putting together the block diagram representing the structure as shown in Figure 5 with the aerodynamic transfer function shown in Figure 7 as a feedback loop based on the structural response to generated aerodynamic forces, the flow structure interaction can be represented as two subsystems interconnected to each other, triggered by a disturbance subsystem as explained in Figure 8.

Figure 8.

Simulation of flow-structure interactions is carried out by arranging the subsystem of structure containing mass, mass moment of inertia, structural stiffness and damping, coupled with aerodynamic forces in lift and pitching moment as functions of the dynamic response of structure.

For validating the numerical model and simulation showed in Figure 8, the case of linear elastic of the aero-elastic system is conducted first. Figures 9 and 10 shows the response of the aerofoil in translation for the airflow speeds of below and above the flutter boundary.

Figure 9.

Convergent translation response of the aerofoil at the flow speed of 13 m/s with initial condition of 4 cm displacement.

Figure 10.

Divergent unstable response of the aerofoil due to a small impulse disturbance at the flow speed of 16 m/s in translation.

The simulation for the aero-elastic system at critical airflow under a certain disturbance or initial condition will evidently yield a constant amplitude of sinusoidal motions.

3.2 Flutter limit cycle oscillation

Phase portraits representing the relationships between the displacement and the velocity of the response are used to analyse the dynamic behaviour of the non-linear system.

Figure 11 presents the phase diagram or the case of linear aero-elastic system where the ellipsoidal trajectories show a stable harmonic response of flutter boundary at the wind speed of flutter boundary.

Figure 11.

Phase portrait of the linear aero-elastic system at the flutter speed boundary of 14 m/s.

A structural non-linear factor may influence the dynamic behaviour of an aero-elastic system. The insertion a free-play in rotation into the support mechanism of the model for an example, at the wind speed below the critical flutter, the system will not reduce entirely the dynamic response of the aerofoil under a certain perturbation as in the linear case, but introduce an oscillation with constant amplitude at a certain frequency. Figure 12 shows a limit cycle oscillation at the wind speed of 13.0 m/s of the aero-elastic system under an initial condition.

Figure 12.

Phase portrait of the aero-elastic system with free-play mechanism of 2° of rotation at the wind speed of 13.0 m/s due to initial displacement of 4 cm.

Furthermore structural non-linear factors may also influence the limit of stability. The free-play mechanism reduces the flutter critical speed for around 0.5 m/s as showed in Figures 13 and 14.

Figure 13.

Phase portrait for the existence of 2° free-play mechanism in rotation at the wind speed of 13.5 m/s with the initial condition of 4 cm displacement.

Figure 14.

Phase portrait for the existence of 2° free-play mechanism in rotation at the wind speed of 14.0 m/s.

At the critical speed calculated for linear system but with the existence of the free-play, the system will generate an unstable divergent structural response as explained in Figure 14.

From the wind-tunnel flutter testing it is observed that a moderate oscillation starts at 13 m/s flow speed and becomes severe vibration at 18 m/s. It can be concluded that a small free-play mechanism involves in the lower speed (less than the critical boundary), and a hardening-stiffness behaviour for the higher speed.

4. Conclusion

A two-degree-of-freedom in transversal and rotational motions wing section for describing classical binary flutter mechanism is used to investigate the effect of free-play nonlinearity to the stability of the aero-elastic system and the associated limit cycles. The aerodynamic forces are calculated by using Theodorsen’s method in frequency domain based on thin aerofoil theory.

By representing the aero-elastic system as a closed loop block diagrams of a dynamic system where the structural part serves as the main plant of the system and the aerodynamic transfer function as a feedback loop calculated based on the dynamic structural response, it is suitable to carry out the simulation on the platform of Simulink-Matlab. For this purpose the aerodynamic forces have to be conversed in Laplace domain.

The work shows the effectiveness of the flow-structure interactions when the system is considered as a dynamic system where the response can be analysed in time domain and the effects of non-linear factors can be conveniently included simultaneously.

The limit cycle oscillation and stability can be showed numerically by representing the phase portrait of the response. At the speed of airflow below the critical speed of flutter, a constant oscillation may happen due to a free-play nonlinearity. It can be shown that the stability boundary becomes smaller than the critical speed.

Acknowledgments

This research is funded by Ministry of Higher Education Malaysia under Fundamental Research Grant Scheme, no. FRGS/2/2013/TK09/02/1.