Hostname: page-component-76c49bb84f-bg8zk Total loading time: 0 Render date: 2025-07-04T02:22:34.476Z Has data issue: false hasContentIssue false

Aeroelastic stability of imperfectly supported high-aspect-ratio wings

Published online by Cambridge University Press:  23 June 2025

M. Kheiri*
Affiliation:
Fluid-Structure Interactions & Aeroelasticity Laboratory, Department of Mechanical, Industrial & Aerospace Engineering, Concordia University, Montreal, QC, Canada
M. Riazat
Affiliation:
Fluid-Structure Interactions & Aeroelasticity Laboratory, Department of Mechanical, Industrial & Aerospace Engineering, Concordia University, Montreal, QC, Canada
*
Corresponding author: M. Kheiri; Email: mojtaba.kheiri@concordia.ca
Rights & Permissions [Opens in a new window]

Abstract

This paper examines the aeroelastic stability of uniform flexible wings imperfectly supported at one end and free at the other. Real-world aircraft wings inevitably exhibit imperfections, including non-ideal end supports. This work is motivated by the critical need to fundamentally understand how end-support imperfections influence the aeroelastic behaviour of fixed wings. The equations of motion are obtained via the extended Hamilton’s principle. The bending-torsional dynamics of the wing is approximated using the Euler-Bernoulli beam theory. The aerodynamic lift and pitching moment are modelled using the unsteady aerodynamics for the arbitrary motion of thin aerofoils in the time domain, extended by the strip flow theory. The imperfect support is modelled via rotational springs (with linear stiffness) for both bending and torsional degrees of freedom. The Galerkin method is used for the spatial discretisation. The stability analysis is performed by solving the resulting eigenvalue problem, and the numerical results are presented in Argand diagrams. The numerical results presented in this study are novel and offer great insights. It is demonstrated that support imperfections can substantially influence the critical flow velocity for both flutter and divergence, as well as alter the sequence of instabilities and the unstable mode. The extent of these effects directly depends on the magnitude of the imperfections. Interestingly–and counterintuitively–in certain cases, a reduction in the flutter speed is observed as the imperfections decrease.

Information

Type
Research Article
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2025. Published by Cambridge University Press on behalf of Royal Aeronautical Society

Nomenclature

Abbreviations

a

dimensionless distance between the elastic axis and mid-chord

b

semi-chord

c

aerofoil chord length

EI

bending stiffness

GJ

torsional stiffness

${I_{ea}}$

mass moment of inertia about the elastic axis

${K_w}$

rotational spring stiffness about the $x$ -axis

${K_\theta }$

rotational spring stiffness about the $y$ -axis

$\ell $

wing span

L

lift

m

mass per unit length

M

aerodynamic pitching moment

t

time

T

kinetic energy

U

flow velocity

V

potential energy

w

displacement in the transverse direction

${x_\theta }$

static unbalance

y

spanwise distance

Greek symbols

$\delta w$

virtual transverse displacement

$\delta \theta $

virtual torsional displacement

$\theta $

pitch angle

$\rho $

fluid density

1.0 Introduction

Aeroelasticity is the science of interactions among aerodynamic, elastic, and inertial loads. These interactions can result into instabilities like divergence and flutter or can manifest themselves, for example, as transient vibrations due to turbulent flow (i.e., buffeting) [Reference Bisplinghoff, Ashley and Halfman1]. Beam-like structures with aerofoil cross-section are found in many engineering systems, such as fixed- and rotary-wing aircraft, wind turbines, compressors and gas turbines. Aeroelastic stability analysis is an essential step in the design process of such systems. For example, according to airworthiness requirements set by the civil aviation authorities, such as the Federal Aviation Administration (FAA), it must be shown that an airplane is free from flutter and divergence for any condition of operation within the flight envelope, including also an adequate (typically 15%) margin of safety [Reference Hodges and Pierce2]. Demonstrating compliance with aeroelastic stability requirements–such as those established by the FAA–typically involves a combination of numerical analyses, wind tunnel testing, ground vibration testing and flight testing [3].

A broad spectrum of aeroelastic models has been developed in the literature, exhibiting substantial variation in aerodynamic and structural dynamic modeling fidelity, as well as in their coupling methodologies. At one end of the spectrum are two- or three-dimensional (3-D) analytical models that utilise aerodynamic theories such as Theodorsen’s thin-aerofoil theory [Reference Bisplinghoff, Ashley and Halfman1] or Peters’ finite-state theory [Reference Peters, Karunamoorthy and Caom4], and represent structural dynamics using a limited number of degrees of freedom. These simplified low-fidelity models are widely used due to their computational efficiency and analytical tractability, making them well-suited for conceptual and exploratory design phases [Reference Afonso, Vale, Oliveira, Lau and Suleman5]. They are additionally very useful for conducting fundamental studies, where understanding the underlying physical mechanisms for the aeroelastic behaviour is the key. Some earlier examples of such aeroelastic models can be found in [Reference Tang and Dowell6Reference Shams, Lahidjani and Haddadpour10] and some more recent examples can be found in [Reference Berci and Cavallaro11Reference Zhang, Kheiri and Xie14].

Many studies have employed 2-D, typical-section aeroelastic models, where ‘distributed’ structural properties such as stiffness are represented using ‘concentrated’ elements. For example, as discussed by Lee et al. [Reference Lee, Price and Wong15], a twisting thin wing or blade behaves like a cubic hardening spring. In other words, a cubic spring in a typical section model can effectively represent the nonlinear stiffness behaviour of a wing undergoing large deformations. Concentrated structural elements such as springs can also be used to represent the structural behaviour in control mechanisms as well as in the parts connecting, for example, a wing to the fuselage or a pylon to the wing. Many studies have focused on the study of the effects of concentrated stiffness nonlinearities on the aeroelastic behaviour (including flutter suppression, dynamic response and control performance) in typical section models; for example, refer to Refs [Reference Zhang, Kheiri and Xie13Reference Amoozgar, Castrichini, Garvey, Friswell, Cooper and Ajaj21]. In fact, the aeroelastic model presented in this paper differs from most studies in the literature by retaining both distributed and concentrated stiffness representations. The distributed stiffness captures the elastodynamic deformation along the wing, while the concentrated stiffness models the structural connection between the wing and the fuselage.

At the other end of the spectrum of aeroelastic models are high-fidelity, strongly (or tightly) coupled computational models that integrate computational fluid dynamics (CFD) and computational structural dynamics (CSD) methods. These models are solved using either a ‘monolithic’ or a ‘partitioned’ scheme. In a monolithic scheme, the CFD and CSD equations are solved simultaneously [Reference Yuan, Sandhu and Poirel22]. Monolithic schemes require reformulation of both fluid and structural dynamics, and they are generally computationally challenging and mathematically suboptimal [Reference Farhat and Lesoinne23], and as a result, they have not received considerable attention. In contrast, ‘partitioned’ schemes allow for spatial and temporal discretisations tailored to the fluid/structure models [Reference Farhat and Lesoinne23]. Despite recent advancements and popularity of computational aeroelastic models (some examples are Refs [Reference Farhat, Geuzaine and Brown24, Reference Dang, Yang and Li25]), they are still time consuming for aeroelastic stability analysis, given that, for example, for an airplane normally hundreds of flight scenarios need to be investigated.

In addition to low- and high-fidelity aeroelastic models, there exist medium-fidelity models, which are commonly used in practice. They typically employ potential-flow-based lifting surface methods, such as the Doublet Lattice Method (DLM) [Reference Albano and Rodden26] or the Vortex Lattice Method (VLM) [Reference Katz and Plotkin27] for aerodynamic modeling. Structural dynamics is generally represented using modal-based methods or finite element models. Among medium-fidelity approaches, DLM-based models are the industry standard for analysing aeroelastic stability and time response in the subsonic flow regime [Reference Afonso, Vale, Oliveira, Lau and Suleman5, Reference Murua, Palacios and Graham28]. They form the foundation of widely used commercial aeroelastic software packages, such as MSC Nastran [29] and ZAERO [30]. Some recent examples of aeroelastic models based on DLM/VLM can be found in [Reference Parenteau and Laurendeau31, Reference Hilger and Ritter32].

Concerning the aeroelasticity of wings, in almost all previous studies, the wing was typically considered to be clamped (or fixed) at one end and free at the other–a cantilevered wing. However, in reality no structure or attachment is perfect. Imperfections may be created, for example, during manufacturing, due to structural fatigue (tear and wear), or during installation, and they may appear in the geometry (e.g., pre-existing deformation, curvature, twist and cross-sectional asymmetry), material (e.g., inhomogeneity, cracks, void mass and delamination) and supports (e.g., loose/flexible end-supports). Studies have indicated the importance of monitoring imperfections and defects to ensure the good health of wing-like structures. For example, blades ‘root attachment problems’ are a common cause of vibration and failure in axial compressors [Reference Meher-Homji and Gabriles33]. In a civil aircraft airframe, fatigue may cause cracks to quickly spread in susceptible structural elements, such as the over-wing fuselage attachment [Reference Munns and Kent34]. This seems particularly crucial as composite materials are becoming increasingly widespread in aerospace and wind energy applications.

Although some research has been conducted in the past to examine the effects of structural damage, such as surface cracks and delamination on the aeroelastic behaviour of wings (e.g. Refs [Reference Wang and Inman35Reference Torabi, Shams, Fatehi Narab and Atashgah38]), the effect of imperfect end-supports on the aeroelastic stability is still unknown. The effects of support imperfection, however, have been examined for some other systems involving fluid-structure interactions, such as cylinders in axial flow [Reference Kheiri and Païdoussis39, Reference Tabatabaei, Kheiri and Dargahi40] and pipes conveying fluid [Reference Kheiri, Païdoussis, Costa Del Pozo and Amabili41Reference Riazat and Kheiri43]. It was found that depending on the system parameters and support imperfection, the critical flow velocity for flutter and the post-flutter dynamical behaviour could significantly be altered. The objective of the present paper is, therefore, to explore the effects of support imperfection on the aeroelastic stability of uniform flexible wings. The rest of the paper is organised as follows. In Section 2.0, the mathematical model is developed. Next, in Section 3.0, the numerical convergence study and the verification of the mathematical model are described. The numerical results are presented and discussed in Section 4.0. And the paper ends with some concluding remarks in Section 5.0.

2.0 Theory

Figure 1. Schematic drawing of a uniform flexible wing imperfectly supported at ${\rm{y}} = 0$ and free at ${\rm{y}} = \ell $ , where $\ell $ is the span. Two rotational springs of the stiffness ${{\rm{K}}_{\rm{w}}}$ and ${{\rm{K}}_{\rm{\theta }}}$ are used to model the support imperfection; the springs are attached to the wing at one end and attached to a rigid support (e.g., fuselage) at the other end; ${\rm{xyz}}$ is the Cartesian coordinate system attached to the undeformed wing; also, ${\rm{c}}$ is the chord length and ${\rm{U}}$ is the freestream velocity.

Figure 1 shows the schematic drawing of a flexible high aspect-ratio wing (i.e. $\ell /c \gg 1$ , $\ell $ and $c$ being the span and chord length of the wing, respectively) subjected to a uniform flow with velocity $U$ . A Cartesian coordinate system is adopted here, where the $x$ -axis is along the chord, pointing from the leading-edge towards the trailing-edge, the $y$ -axis is along the span and the $z$ -axis along the wing thickness. The wing has uniform mass per unit length $m$ , spanwise bending rigidity $EI$ and torsional rigidity $GJ$ . It is assumed that the $y$ -axis is coincident with the elastic axis of the wing. The bending-torsional dynamics of the wing is modelled using the Euler-Bernoulli beam theory, where $w\left( {y,t} \right)$ denotes bending displacements in the $z$ -direction while $\theta \left( {y,t} \right)$ represents torsional displacements in the $y$ -direction; $t$ represents time.

The wing is imperfectly supported at $y = 0$ while it is free at $y = \ell $ . The support imperfections are considered for both bending and torsional degrees-of-freedom (DOFs), using rotational springs of linear stiffness. Modelling support imperfections using linear springs is a simplification, as imperfections in real-world airframe structures can introduce both stiffness nonlinearities and dissipative effects [Reference Lee, Price and Wong15]. Nevertheless, the present formulation is the first theoretical modelling attempt, and it can easily be extended to account for these complexities. Examples of such extensions for a pipe conveying fluid system are the works by Kheiri [Reference Kheiri42] and Riazat and Kheiri [Reference Riazat and Kheiri43] in which cubic nonlinearities were used to model support imperfections.

The stiffness of the springs is represented by ${K_w}$ and ${K_\theta }$ , respectively. In other words, at $y = 0$ , $w = 0$ while $\partial w/\partial y \ne 0$ and $\theta \ne 0$ . Following the approach introduced by Kheiri et al. [Reference Kheiri, Païdoussis, Costa Del Pozo and Amabili41], the effect of support imperfections is considered in the equations of motion rather than the boundary conditions.

The equations of motion are obtained via the extended Hamilton’s principle, following the formulation presented in Ref. [Reference Hodges and Pierce2]:

(1) \begin{align}T = \frac{1}{2}\mathop \smallint \nolimits_0^\ell \left( {m{{\dot w}^2} - 2mb{x_\theta }\dot w\dot \theta + {I_{ea}}{{\dot \theta }^2}} \right)\mathrm{dy},\end{align}

where $T$ is the kinetic energy of the wing; $b$ is the semi-chord (i.e. $c = 2b$ ), ${x_\theta } = e - a$ is the so-called static unbalance parameter, and ${I_{ea}}$ is the mass moment of inertia per unit length about the elastic axis; also, (.) $ = \partial \left( {{\rm{\;\;}}} \right)/\partial t$ .

The potential energy may be written as

(2) \begin{align}V = \frac{1}{2}\mathop \smallint \nolimits_0^\ell \left( {EI{{w''}^2} + GJ{{\theta '}^2}} \right)\mathrm{dy} + \frac{1}{2}{K_w}{w'^2}\left( {0,t} \right) + \frac{1}{2}{K_\theta }{\theta ^2}\left( {0,t} \right),\end{align}

where the integral term represents the strain energy due to bending and torsion and the last two terms on the r.h.s. denote the potential energy due to springs (used to model support imperfections); also, ${({\rm{\;\;}})^{\prime}} = \partial \left( {{\rm{\;\;}}} \right)/\partial y$ .Footnote 1

Thus, the Lagrangian of the system, i.e., ${\mathcal L} = T - V$ , is obtained and its variation is added to the virtual work due to aerodynamic lift per unit length $L$ and pitching moment (about elastic axis) per unit length ${M_{ea}}$ :

(3) \begin{align}\delta {W_{nc}} = \mathop \smallint \nolimits_0^\ell \left( {L\delta w + {M_{ea}}\delta \theta } \right)\mathrm{dy},\end{align}

where $\delta w$ and $\delta \theta $ are the bending and torsional virtual displacements, respectively.

After some mathematical manipulations, the equations of motion in the differential form can be written as

(4) \begin{align} {}& \mathop \int \nolimits_0^\ell \left\{ {\big[ {m\ddot w - mb{x_\theta }\ddot \theta + EI{w^{{\rm{\prime\prime\prime\prime}}}} - L\big] {\delta w} + \big[{K_w}w'\bar \delta ( y)} \big]\delta w'} \right\}\mathrm{dy} = 0,\nonumber\\ {}&\mathop \int \nolimits_0^\ell \big[ {{I_{ea}}\ddot \theta - mb{x_\theta }\ddot w - GJ\theta '' + {K_\theta }\theta \bar \delta ( y) - {M_{ea}}} \big]\delta \theta\, \mathrm{dy} = 0,\end{align}

where $\bar \delta \left( y \right)$ denotes the Dirac delta function.

In this paper, the indicial aerodynamic theory based on Wagner’s function is used to represent the unsteady aerodynamic lift and pitching moment in the time domain for small arbitrary motions of the wing. The two-dimensional aerodynamic theory is extended to three dimensions using the strip flow theory. This approach is widely used in the literature; some recent examples are Refs [Reference Asadi and Farsadi44Reference Vindigni, Mantegna, Orlando and Alaimo47]. Since the indicial aerodynamic theory is valid for arbitrary motions, as opposed to Theodorsen’s thin-aerofoil theory, which is valid for simple harmonic motions, the solution is accurate not only at the onset of flutter, but also for all flow velocities below the flutter speed. Additionally, a time domain formulation facilitated by the use of the indicial aerodynamic theory can be used for the investigation of the dynamic response of the aeroelastic system to external disturbances such as gusts. Moreover, as seen in the following, the stability analysis is conducted more straightforwardly, compared to an iterative solution process, which is typical of frequency domain stability analyses.

It is noted that since the wing is assumed to be of high-aspect ratio, the effects of three-dimensionality of the flow and tip vortices are neglected. Ignoring such effects typically leads to aerodynamic loads that are higher than those encountered in reality, which in turn results in a lower flutter speed – a more conservative prediction. A recent comparative study [Reference Modaress-Aval, Bakhtiari-Nejad, Dowell, Peters and Shahverdi48] showed that for wings with the aspect ratio of 15 and higher, a two-dimensional aerodynamic model is sufficiently accurate to conduct an aeroelastic analysis.

According to Wagner’s problem of the step change in angle-of-attack, lift may be written as [Reference Bisplinghoff, Ashley and Halfman1]:

(5) \begin{align}L &= \pi \rho {b^2}\left[ { - \ddot w + U\dot \theta - ba\ddot \theta } \right] \nonumber\\ &\quad- 2\pi \rho Ub\left( {{v_{3/4}}\left( 0 \right)\varphi \left( \tau \right) + \mathop \smallint \nolimits_0^\tau \frac{{{\mathrm{d}}{v_{3/4}}\left( \sigma \right)}}{{{\mathrm{d}}\sigma }}\varphi \left( {\tau - \sigma } \right){\mathrm{d}}\sigma } \right),\end{align}

where ${v_{3/4}}\left( t \right) = \dot w - U\theta - b\left( {1/2 - a} \right)\dot \theta $ is the instantaneous vertical velocity of the fluid particle (the so-called downwash) in contact with the three-quarter chord point of the aerofoil section; $a$ is the dimensionless distance between the elastic axis and mid-chord; $\tau = Ut/b$ is the dimensionless time; also, $\varphi \left( \tau \right)$ represents Wagner’s function, which can conveniently be approximated as

(6) \begin{align}\varphi \left( \tau \right) = 1 - {\gamma _1}{e^{ - {\varepsilon _1}\tau }} - {\gamma _2}{e^{ - {\varepsilon _2}\tau }},\end{align}

in which ${\gamma _1} = 0.165$ , ${\gamma _2} = 0.335$ , ${\varepsilon _1} = 0.0455$ , and ${\varepsilon _2} = 0.3$ .

The aerodynamic pitching moment about the elastic axis may also be written as

(7) \begin{align}{M_{ea}} &= \pi \rho {b^2}\left[ { - ba\ddot w - Ub\left( {\frac{1}{2} - a} \right)\dot \theta - {b^2}\left( {\frac{1}{8} + {a^2}} \right)\ddot \theta } \right]\nonumber \\ &\quad- 2\pi \rho U{b^2}\left( {\frac{1}{2} + a} \right)\left( {{v_{3/4}}\left( 0 \right)\varphi \left( \tau \right) + \mathop \smallint \nolimits_0^\tau \frac{{{\mathrm{d}}{v_{3/4}}\left( \sigma \right)}}{{{\mathrm{d}}\sigma }}\varphi \left( {\tau - \sigma } \right){\mathrm{d}}\sigma } \right).\end{align}

By substituting Equation (6) into the expression for lift (see Equation (5)), we obtain:

(8) \begin{align} L &= \pi \rho {b^2}\left[ { - \ddot w + U\dot \theta - ba\ddot \theta \left] { - 2\pi \rho Ub} \right[\dot w - U\theta - b\left( {\frac{1}{2} - a} \right)\dot \theta } \right]\varphi \left( 0 \right)\nonumber\\ &\quad- 2\pi \rho {U^2}\left[ {w - b\left( {\frac{1}{2} - a} \right)\theta } \right]\frac{{{\mathrm{d}}\varphi }}{{{\mathrm{d}}\tau }}\left( 0 \right) + 2\pi \rho {U^2}\left[ {w\left( 0 \right) - b\left( {\frac{1}{2} - a} \right)\theta \left( 0 \right)} \right]\frac{{{\mathrm{d}}\varphi }}{{{\mathrm{d}}\tau }}\nonumber \\ &\quad+ 2\pi \rho {U^2}\mathop \sum \limits_{k = 1}^2 \left[ {b{\gamma _k}{\varepsilon _k}\left( {1 - \left( {\frac{1}{2} - a} \right){\varepsilon _k}} \right){A_k} + {\gamma _k}\varepsilon _k^2{H_k}} \right],\end{align}

and the expression for the pitching moment becomes

(9) \begin{align} {M_{ea}} &= \pi \rho {b^2}\left[ { - ba\ddot w - Ub\left( {\frac{1}{2} - a} \right)\dot \theta - {b^2}\left( {\frac{1}{8} + {a^2}} \right)\ddot \theta } \right]\nonumber\\&\quad - 2\pi \rho U{b^2}\left( {\frac{1}{2} + a} \right)\left[ {\dot w - U\theta - b\left( {\frac{1}{2} - a} \right)\dot \theta } \right]\varphi \left( 0 \right)\nonumber\\&\quad - 2\pi \rho {U^2}b\left( {\frac{1}{2} + a} \right)\left[ {w - b\left( {\frac{1}{2} - a} \right)\theta } \right]\frac{{d\varphi }}{{d\tau }}\left( 0 \right)\nonumber\\&\quad + 2\pi \rho {U^2}b\left( {\frac{1}{2} + a} \right)\left[ {w\left( 0 \right) - b\left( {\frac{1}{2} - a} \right)\theta \left( 0 \right)} \right]\frac{{d\varphi }}{{d\tau }} \nonumber \\&\quad + 2\pi \rho {U^2}b\left( {\frac{1}{2} + a} \right)\mathop \sum \limits_{k = 1}^2 \left[ {b{\gamma _k}{\varepsilon _k}\left( {1 - \left( {\frac{1}{2} - a} \right){\varepsilon _k}} \right){A_k} + {\gamma _k}\varepsilon _k^2{H_k}} \right],\end{align}

where ${A_k}$ and ${H_k}$ ( $k = 1,2)$ are:

(10) \begin{align}{A_k} = \mathop \smallint \nolimits_0^\tau {e^{ - {\varepsilon _k}\left( {\tau - \sigma } \right)}}\theta \left( \sigma \right){\mathrm{d}}\sigma ,{H_k} = \mathop \smallint \nolimits_0^\tau {e^{ - {\varepsilon _k}\left( {\tau - \sigma } \right)}}w\left( \sigma \right){\mathrm{d}}\sigma .\end{align}

Following dimensionless variables are defined

(11) \begin{align}& \eta = \frac{w}{\ell },{\rm{\;\;\;\;}}\zeta = \frac{y}{\ell },{\rm{\;\;\;\;}}\tau = \frac{{Ut}}{b},{\rm{\;\;\;\;}}\mu = \frac{m}{{\pi \rho {b^2}}},{\rm{\;\;\;\;}}u = {\left(\frac{{\pi \rho {b^2}}}{{EI}}\right)^{\frac{1}{2}}}U\ell ,\nonumber\\& AR = \frac{\ell }{b},{\rm{\;\;\Gamma }} = \frac{{EI}}{{GJ}},{\rm{\;\;\;\;}}{r^2} = \frac{{{I_{ea}}}}{{m{b^2}}},{\rm{\;\;\;\;}}{k_\eta } = \frac{{{K_w}{\ell ^3}}}{{EI}},{\rm{\;\;\;\;}}{k_\theta } = \frac{{{K_\theta }\ell }}{{GJ}},\end{align}

where $\mu $ is called mass ratio, $AR$ aspect ratio, ${\rm{\Gamma }}$ stiffness ratio and $r$ the dimensionless radius of gyration.

Using the above dimensionless variables, the partial differential equations (Equation (4)) are rendered dimensionless. Then, Galerkin’s method is utilised to spatially discretise the dimensionless partial differential equations by assuming $\eta \left( {\zeta ,\tau } \right) = {\rm{\Sigma }}_{i = 1}^N{{\rm{\Phi }}_i}\left( \zeta \right){q_i}\left( \tau \right)$ and $\theta \left( {\zeta ,\tau } \right) = {\rm{\Sigma }}_{i = 1}^N{{\rm{\Theta }}_i}\left( \zeta \right){p_i}\left( \tau \right)$ , where ${{\rm{\Phi }}_i}\left( \zeta \right)$ and ${{\rm{\Theta }}_i}\left( \zeta \right)$ are, respectively, the bending and torsional mode shapes, and ${q_i}\left( \tau \right)$ and ${p_i}\left( \tau \right)$ are their corresponding generalised coordinates; also, $N$ is the number of modes, which is assumed to be the same for bending and torsion.

The discretised, dimensionless bending-torsional equations of motion are expressed as

(12) \begin{align} &M_{ij}^{\left( 1 \right)}{\ddot q_i} + M_{ij}^{\left( 2 \right)}{\ddot p_i} + C_{ij}^{\left( 1 \right)}{\dot q_i} + C_{ij}^{\left( 2 \right)}{\dot p_i} + K_{ij}^{\left( 1 \right)}{q_i} + K_{ij}^{\left( 2 \right)}{p_i} + G_{ij}^{\left( 1 \right)}{H_{1i}} + G_{ij}^{\left( 2 \right)}{H_{2i}} \nonumber\\ &\quad+ G_{ij}^{\left( 3 \right)}{A_{1i}} + G_{ij}^{\left( 4 \right)}{A_{2i}} + F_i^{\left( 1 \right)}\left( \tau \right) = 0,\end{align}

and

(13) \begin{align} &M_{ij}^{\left( 3 \right)}{\ddot q_i} + M_{ij}^{\left( 4 \right)}{\ddot p_i} + C_{ij}^{\left( 3 \right)}{\dot q_i} + C_{ij}^{\left( 4 \right)}{\dot p_i} + K_{ij}^{\left( 3 \right)}{q_i} + K_i^{\left( 4 \right)}{p_i} + G_{ij}^{\left( 5 \right)}{H_{1i}} + G_{ij}^{\left( 6 \right)}{H_{2i}}\nonumber\\ &\quad+ G_{ij}^{\left( 7 \right)}{A_{1i}} + G_{ij}^{\left( 8 \right)}{A_{2i}} + F_i^{\left( 2 \right)}\left( \tau \right) = 0,\end{align}

where

(14) \begin{align}M_{ij}^{\left( 1 \right)} &= \left( {1 + \mu } \right){u^2}A{R^2}{\delta _{ij}},{\rm{\;\;}}M_{ij}^{\left( 2 \right)} = \left( {a - \mu {x_\theta }} \right){u^2}AR{b_{ij}},{\rm{\;\;}}C_{ij}^{\left( 1 \right)} = 2{u^2}A{R^2}\varphi \left( 0 \right){\delta _{ij}},\nonumber\\C_{ij}^{\left( 2 \right)} &= - {u^2}AR\left[ {1 + 2\left( {\frac{1}{2} - a} \right)\varphi \left( 0 \right)} \right]{b_{ij}},{\rm{\;\;}}K_{ij}^{\left( 1 \right)} = \lambda _i^4{\delta _{ij}} + 2{u^2}A{R^2}\dot \varphi \left( 0 \right){\delta _{ij}} + {k_\eta }{\rm{\Phi }}_i^{\rm{\prime}}\left( 0 \right){\rm{\Phi }}_j^{\rm{\prime}}\left( 0 \right)\nonumber\\K_{ij}^{\left( 2 \right)} &= - 2{u^2}AR\left[ {\varphi \left( 0 \right) + \left( {\frac{1}{2} - a} \right)\dot \varphi \left( 0 \right)} \right]{b_{ij}}\nonumber\\G_{ij}^{\left( 1 \right)} &= - 2{u^2}A{R^2}{\gamma _1}\varepsilon _1^2{\delta _{ij}},{\rm{\;\;}}G_{ij}^{\left( 2 \right)} = - 2{u^2}A{R^2}{\gamma _2}\varepsilon _2^2{\delta _{ij}},\nonumber\\G_{ij}^{\left( 3 \right)} &= - 2{u^2}AR{\gamma _1}{\varepsilon _1}\left( {1 - \left( {\frac{1}{2} - a} \right){\varepsilon _1}} \right){b_{ij}},{\rm{\;\;}}G_{ij}^{\left( 4 \right)} = - 2{u^2}AR{\gamma _2}{\varepsilon _2}\left( {1 - \left( {\frac{1}{2} - a} \right){\varepsilon _2}} \right){b_{ij}}\nonumber\\F_i^{\left( 1 \right)} &= - 2{u^2}A{R^2}{\delta _{ij}}{q_i}\left( 0 \right)\dot \varphi \left( \tau \right) + 2{u^2}AR\left( {\frac{1}{2} - a} \right){b_{ij}}{p_i}\left( 0 \right)\dot \varphi \left( \tau \right),\end{align}
(15) \begin{align}M_{ij}^{\left( 3 \right)} &= \left( {a - \mu {x_\theta }} \right){\rm{\Gamma }}{u^2}AR{b_{ji}},{\rm{\;\;}}M_{ij}^{\left( 4 \right)} = \left( {\mu {r^2} + \frac{1}{8} + {a^2}} \right){\rm{\Gamma }}{u^2}{\delta _{ij}},{\rm{\;\;}}C_{ij}^{\left( 3 \right)} = 2{\rm{\Gamma }}{u^2}AR\left( {\frac{1}{2} + a} \right)\varphi \left( 0 \right){b_{ji}},\nonumber\\C_{ij}^{\left( 4 \right)} &= {\rm{\Gamma }}{u^2}\left( {\frac{1}{2} - a} \right)\left[ {1 - 2\left( {\frac{1}{2} + a} \right)\varphi \left( 0 \right)} \right]{\delta _{ij}},{\rm{\;\;}}K_{ij}^{\left( 3 \right)} = 2{\rm{\Gamma }}{u^2}AR\left( {\frac{1}{2} + a} \right)\dot \varphi \left( 0 \right){b_{ji}},\nonumber\\K_{ij}^{\left( 4 \right)} &= - \left[ {{c_{ij}} + 2{\rm{\Gamma }}{u^2}\left( {\frac{1}{2} + a} \right)\varphi \left( 0 \right){\delta _{ij}} + 2{\rm{\Gamma }}{u^2}\left( {\frac{1}{4} - {a^2}} \right)\dot \varphi \left( 0 \right){\delta _{ij}} - {k_\theta }{{\rm{\Theta }}_i}\left( 0 \right){{\rm{\Theta }}_j}\left( 0 \right)} \right]\nonumber\\G_{ij}^{\left( 5 \right)} &= - 2{\rm{\Gamma }}{u^2}AR\left( {\frac{1}{2} + a} \right){\gamma _1}\varepsilon _1^2{b_{ji}},{\rm{\;\;}}G_{ij}^{\left( 6 \right)} = - 2{\rm{\Gamma }}{u^2}AR\left( {\frac{1}{2} + a} \right){\gamma _2}\varepsilon _2^2{b_{ji}},\nonumber\\G_{ij}^{\left( 7 \right)} &= - 2{\rm{\Gamma }}{u^2}\left( {\frac{1}{2} + a} \right){\gamma _1}{\varepsilon _1}\left( {1 - \left( {\frac{1}{2} - a} \right){\varepsilon _1}} \right){\delta _{ij}},{\rm{\;\;}}G_{ij}^{\left( 8 \right)} = - 2{\rm{\Gamma }}{u^2}\left( {\frac{1}{2} + a} \right){\gamma _2}{\varepsilon _2}\left( {1 - \left( {\frac{1}{2} - a} \right)\!{\varepsilon _2}} \right)\!{\delta _{ij}}\nonumber\\F_i^{\left( 2 \right)} &= - 2{\rm{\Gamma }}{u^2}AR\left( {\frac{1}{2} + a} \right){b_{ji}}{q_i}\left( 0 \right)\dot \varphi \left( \tau \right) + 2{\rm{\Gamma }}{u^2}\left( {\frac{1}{4} - {a^2}} \right){\delta _{ij}}{p_i}\left( 0 \right)\dot \varphi \left( \tau \right),\end{align}

in which ${\delta _{ij}} = \mathop \smallint \nolimits_0^1 {{\rm{\Phi }}_i}{{\rm{\Phi }}_j}{\mathrm{d}}\zeta = \mathop \smallint \nolimits_0^1 {{\rm{\Theta }}_i}{{\rm{\Theta }}_j}{\mathrm{d}}\zeta $ is the Kronecker’s delta function (considering orthonormal mode shapes ${{\rm{\Phi }}_i}$ and ${{\rm{\Theta }}_i}$ ), ${b_{ij}} = \mathop \smallint \nolimits_0^1 {{\rm{\Theta }}_i}{{\rm{\Phi }}_j}{\mathrm{d}}\zeta $ , and ${c_{ij}} = \mathop \smallint \nolimits_0^1 {\rm{\Theta }}_i^{{\rm{\prime\prime}}}{{\rm{\Theta }}_j}{\mathrm{d}}\zeta $ ; ${A_{ki}} = \mathop \smallint \nolimits_0^\tau {e^{ - {\varepsilon _k}\left( {\tau - \sigma } \right)}}{p_i}{\mathrm{d}}\sigma $ and ${H_{ki}} = \mathop \smallint \nolimits_0^\tau {e^{ - {\varepsilon _k}\left( {\tau - \sigma } \right)}}{q_i}{\mathrm{d}}\sigma $ ; note that the fact that ${\rm{\Phi }}_i^{{\rm{\prime\prime\prime\prime}}} = \lambda _i^4{{\rm{\Phi }}_i}$ has been employed (see below); also, the overdot and prime have been re-defined as (.) $ = d\left( {{\rm{\;\;}}} \right)/d\tau $ and ${({\rm{\;\;}})^{\rm{\prime}}} = d\left( {{\rm{\;\;}}} \right)/d\zeta $ .

So far in the modelling, imperfections have been considered for both bending and torsional DOFs; however, for simplicity, in the parametric study and in the analysis, we consider imperfections only for the torsional DOF, i.e. ${k_\eta } \to \infty $ . Therefore, ${{\rm{\Phi }}_i}(\zeta $ ) represent mode shapes for free vibration of a clamped-free beam in bending, while ${{\rm{\Theta }}_i}\left( \zeta \right)$ represent the mode shapes for free vibration of a free-free beam in torsion:

(16) \begin{align}&{{\rm{\Phi }}_i}\left( \zeta \right) = {\mathrm{cosh}}{\lambda _i}\zeta - {\mathrm{cos}}{\lambda _i}\zeta - {\beta _i}\left( {{\mathrm{sinh}}{\lambda _i}\zeta - {\mathrm{sin}}{\lambda _i}\zeta } \right),{\rm{\;\;\;\;}}i = 1 \cdots N\end{align}
(17) \begin{align}&{{\rm{\Theta }}_1}\left( \zeta \right) = 1,\nonumber\\&{{\rm{\Theta }}_i}\left( \zeta \right) = \sqrt 2 {\mathrm{cos}}{\alpha _i}\zeta ,{\rm{\;\;\;\;}}i = 2 \cdots N\end{align}

where ${\lambda _i}$ , ${\beta _i}$ and ${\alpha _i}$ are given in Ref. [Reference Hodges and Pierce2]; also, ${{\rm{\Theta }}_1}$ represents the rigid-body mode shape. Note that ${{\rm{\Phi }}_i}$ and ${{\rm{\Theta }}_i}$ are orthonormal, i.e. $\mathop \smallint \nolimits_0^1 {{\rm{\Phi }}_i}{{\rm{\Phi }}_i}{\mathrm{d}}\zeta = \mathop \smallint \nolimits_0^1 {{\rm{\Theta }}_i}{{\rm{\Theta }}_i}{\mathrm{d}}\zeta = 1$ and $\mathop \smallint \nolimits_0^1 {{\rm{\Phi }}_i}{{\rm{\Phi }}_j}{\mathrm{d}}\zeta = \mathop \smallint \nolimits_0^1 {{\rm{\Theta }}_i}{{\rm{\Theta }}_j}{\mathrm{d}}\zeta = 0,{\rm{\;\;}}\left( {i \ne j} \right)$ .

Equations (12) and (13) can be cast into the first-order or state-space form by defining a new vector of aeroelastic states ${\mathbf{X}} = {[{\dot{\mathbf{q}}}{\rm{\;\;}}{\dot{\mathbf{p}}}{\rm{\;\;}}{\mathbf{q}}{\rm{\;\;}}{\mathbf{p}}{\rm{\;\;}}{{\mathbf{H}}_1}{\rm{\;\;}}{{\mathbf{H}}_2}{\rm{\;\;}}{{\mathbf{A}}_1}{\rm{\;\;}}{{\mathbf{A}}_2}]^T}$ , where, for example, ${\dot{\mathbf{q}}} = {\{ {\dot q_1}{\dot q_2} \cdots {\dot q_N}\} ^T}$ , and ${{\mathbf{H}}_1} = {\{ {H_{11}}{H_{12}} \cdots {H_{1N}}\} ^T}$ :

(18) \begin{align}{\mathcal A}{\dot{\mathbf{X}}} = {\mathcal B}{\mathbf{X}} + {\mathbf{F}}\left( \tau \right),\end{align}

in which matrices ${\mathcal A}$ , ${\mathcal B}$ and vector ${\mathbf{F}}\left( \tau \right)$ are provided in Appendix 6. It should also be noted that using the Leibniz integral rule, one can show that ${{\dot{\mathbf{H}}}_{\mathbf{k}}} = {\mathbf{q}} - {\varepsilon _k}{{\mathbf{H}}_{\mathbf{k}}}$ and ${{\dot{\mathbf{A}}}_{\mathbf{k}}} = {\mathbf{p}} - {\varepsilon _k}{{\bf{A}}_{\mathbf{k}}}$ .

Assuming ${q_i}\left( 0 \right) = {p_i}\left( 0 \right) = 0$ and ${\mathbf{X}} = {\bar {\mathbf {X}}}{\mathrm{exp}}i\omega \tau $ , where $\omega $ called eigenfrequency, Equation (18) is re-formulated as an eigenvalue problem which can then be solved for stability analysis, using, for example, eig function of MATLAB.

3.0 Convergence Study and Verification of the Model

In order to find the minimum number of mode shapes required for obtaining accurate numerical results, we obtained the dimensionless flow velocity for flutter and divergence, represented by ${u_{cf}}$ and ${u_{cd}}$ , respectively, for different numbers of mode shapes used in the Galerkin approximation for an imperfectly supported wing with ${k_\theta } = {10^6}$ – essentially a cantilevered wing. The parameters for the wing of a high-altitude long-endurance (HALE) aircraft [Reference Patil, Hodges and Cesnik7] were used, which are given in Table 1.

Table 1. Parameters of the wing used in Ref. [Reference Patil, Hodges and Cesnik7] and the corresponding dimensionless parameters

As seen from Table 2, with 10 modes for bending and the same number of modes for torsion, the results are within an acceptable range ( $ \le 1{\rm{\% }}$ ) with respect to those with 14 modes. These numbers of modes were also found to yield accurate numerical results for imperfectly supported wings with finite ${k_\theta }$ ; for example, for a wing with ${k_\theta } = 1$ , even with $N = 6$ , the relative errors for ${u_{cd}}$ and ${u_{cf}}$ were below 1%. However, to ensure greater accuracy, $N = 14$ is adopted for the rest of numerical solutions in this paper.

Table 2. Numerical convergence study using a different number of mode shapes for a wing with the parameters given in Table 1; also, ${{\rm{k}}_{\rm{\theta }}} = {10^6}$

To verify (following the definition provided in [49] for verification) the present model, the critical speeds for the aeroelastic instabilities of two different wings are numerically obtained and compared with those reported in the literature. The first wing is the wing of a HALE [Reference Patil, Hodges and Cesnik7] (see Table 1), and the other is the Goland wing [Reference Goland50] whose parameters can also be found in [Reference Haddadpour and Firouz-Abadi51]. As seen from Table 3, the numerical results from the present model are in very good agreement with those reported in the literature; the maximum relative error is below 5%.

Table 3. Comparison between the critical speeds for the aeroelastic instabilities of two wings obtained from the present model and those found from the literature

4.0 Results and Discussion

As discussed previously, in the present numerical studies, the support imperfection is only considered for the torsional DOF at $\zeta = 0$ , and the wing is considered perfectly clamped at $\zeta = 0$ as far as bending vibrations are concerned. This is done to isolate the impact of support imperfections on the stability of the system, while reducing the influence of interactions between different imperfections.

Here, the dynamics of imperfectly supported wings are described by showing the evolution of the first few eigenfrequencies of the system, represented by $\omega $ , in Argand diagrams, as the dimensionless flow velocity $u$ is varied. In the Argand diagram, which is a complex plane, the abscissa and ordinate correspond to the real and imaginary parts of $\omega $ ( $Re\left( \omega \right)$ and $Im\left( \omega \right)$ ), respectively; $Re\left( \omega \right)$ corresponds to the dimensionless frequency while $Im\left( \omega \right)$ is related to damping. In the Argand diagram, the positive half-plane ( $Im\left( \omega \right) \gt 0$ ) is the stable half-plane whereas the negative half-plane ( $Im\left( \omega \right) \lt 0$ ) is the unstable half-plane. If an eigenfrequency locus crosses from the positive half-plane to the negative half-plane while $Re\left( \omega \right) \gt 0$ , the instability is flutter, whereas if $Re\left( \omega \right) = 0$ , the instability is divergence; for more details, the reader is referred to Ref. [Reference Païdoussis52].

Figure 2(a–d) shows the Argand diagrams for a wing with the parameters of the HALE wing (Table 1) as the imperfection becomes more significant (i.e. ${k_\theta }$ is reduced). In the plots, the numerals printed along each locus represent values of $u$ . Figure 2(a) shows the Argand diagram for the wing with almost no imperfection ( ${k_\theta } = {10^6}$ ). As seen, for $u \gt 0$ , all modes (or eigenfrequencies) lie in the stable half-plane, which means that flow induces positive damping in all modes and any vibrations due to external disturbances would die out with time. However, as $u$ is increased, damping in the first torsional mode gradually diminishes, and the locus eventually crosses from the stable half-plane to the unstable half-plane at $u = {u_{cf}} \simeq 1.01$ , while $Re\left( \omega \right) \ne 0$ , indicating flutter. The mode evolving on the $Im\left( \omega \right)$ -axis also becomes unstable at a slightly higher flow velocity, i.e. $u = {u_{cd}} \simeq 1.13$ , indicating divergence. This mode, which is called ‘zeroth mode’ in the plots, is due to either ${A_k}$ or ${H_k}$ generalised coordinates. In fact, 4 $N$ loci with zero frequency are obtained numerically; however, to have clear Argand diagrams, the stable loci have been filtered out and only the unstable loci have been shown. From the plot, it is seen that the first torsional mode is re-stabilised at $u \simeq 2.40$ while the system remains statically unstable due to the zeroth mode. The other modes remain stable up to the maximum flow velocity investigated, i.e. ${u_{max}} = 3$ .

Figure 2. Argand diagrams showing the evolution of the first few dynamical modes of an imperfectly-supported wing with parameters given in Table 1 for: (a) ${{\rm{k}}_{\rm{\theta }}} = {10^6}$ , (b) ${{\rm{k}}_{\rm{\theta }}} = {10^0}$ , (c) ${{\rm{k}}_{\rm{\theta }}} = {10^{ - 1}}$ , and (d) ${{\rm{k}}_{\rm{\theta }}} = {10^{ - 6}}$ . The numerals in the plots, close to the loci correspond to the values of the dimensionless flow velocity. The ${\rm{x}}$ and ${\rm{y}}$ axes have been normalised with respect to the natural frequency of the first torsional mode of a cantilevered beam; also, ${{\rm{u}}_{{\rm{cd}}}}$ and ${{\rm{u}}_{{\rm{cf}}}}$ represent the critical value of ${\rm{u}}$ for divergence and flutter, respectively.

It should be emphasised that the modes or eigenfrequencies shown in the Argand diagrams are aeroelastic modes rather than classical vibration modes of a beam. This means that each locus contains appreciable contents from various vibration modes and the aerodynamics. Nevertheless, for ease of discussion, here, we have adopted a similar appellation as that used in Ref. [Reference Païdoussis52] for the loci of eigenfrequencies. A locus of eigenfrequency is called the ‘1st bending mode’ since its frequency (i.e., $Re\left( \omega \right)$ ) at $u = 0$ matches that of the beam’s uncoupled, first bending-mode frequency, and a locus is called the ‘1st torsional mode’ because $Re\left( \omega \right)$ at $u = 0$ matches the frequency of the first torsional mode of the beam, etc.

As seen from Fig. 2(b), in contrast to the wing with almost no support imperfection, for ${k_\theta } = {10^0}$ (i.e. a wing with significant support imperfections), the stability is lost by divergence at a much lower flow velocity ( $u = {u_{cd}} \simeq 0.62$ ) and by flutter at $u = {u_{cf}} \simeq 1.98$ . While the system is already unstable, the second bending mode also undergoes flutter at $u \simeq 2.37$ . Note that by decreasing the end-spring stiffness (used to model support imperfection), the frequency of the torsional modes is generally diminished.

By further increasing the support imperfection (i.e. by decreasing ${k_\theta }$ to ${k_\theta } = {10^{ - 1}}$ ), as shown in Fig. 2(c), the first torsional mode of the system loses stability by flutter at $u = {u_{cf}} \simeq 0.22$ while the mode evolving on the $Im\left( \omega \right)$ -axis becomes unstable at a slightly higher flow ( $u = {u_{cd}} \simeq 0.23$ ). Higher modes, such as second torsional and bending modes, also become unstable at higher flows.

Finally, Fig. 2(d) shows the Argand diagram for the system with effectively no torsional stiffness at $\zeta = 0$ . As seen from the plot, the system loses stability by divergence at essentially $u = {0^ + }$ . At higher values of $u$ , the system loses stability by flutter in the second torsional mode at $u = {u_{cf}} \simeq 1.77$ and in the second bending mode at $u \simeq 2.42$ .

As seen from the Argand diagrams, by increasing imperfections, the dynamical behaviour may significantly change. This includes the reduction in the critical flow velocities for the instabilities, the change in the sequence of the instabilities (e.g. first flutter and then divergence or first divergence and then flutter), as well as the unstable mode. Nevertheless, it must be noted that all of the above predictions are based on the linear theory developed in Section 2.0. It is well-known that a linear theory remains accurate until the first instability is reached, and past the threshold of the first instability, predictions made by linear theories may not be trusted. A nonlinear theory should be developed in the future to verify the predictions made by the linear theory.

Next, we examine the sensitivity of the flutter speed of the imperfectly-supported wing to the system parameters, such as mass ratio, stiffness ratio and the dimensionless radius of gyration. Figure 3 shows the variation of ${u_{cf}}$ as a function of ${k_\theta }$ for various mass ratios. Additionally, the variation of the dimensionless divergence speed (dashed line), which is independent of $\mu $ , is shown. As seen from the plot, when ${k_\theta } \to \infty $ or ${k_\theta } \to 0$ , by increasing $\mu $ , ${u_{cf}}$ decreases. This might look counterintuitive at first glance since one generally expects to see the flutter speed increasing with the increase of mass. This is, in fact, confirmed if one obtains the dimensional flutter speed ${U_f}$ ; thus, one should regard the difference between the dimensionless and dimensional parameters. Also, note that the dimensionless flow velocity was defined differently in this paper compared to what is generally adopted in the literature for typical section models.

Figure 3. Variation of the dimensionless critical flow velocity for flutter ( ${{\rm{u}}_{{\rm{cf}}}}$ ) of an imperfectly supported wing as a function of the end-spring stiffness ( ${{\rm{k}}_{\rm{\theta }}}$ ) for different values of mass ratio. Also, the dashed line shows the critical flow velocity for divergence, which is independent of ${\rm{\mu }}$ . The rest of system parameters are the same as the dimensionless variables given in Table 1.

As also seen from Fig. 3, for each $\mu $ , ${u_{cf}}$ for ${k_\theta } \to 0$ (i.e. pinned-free wing) is larger than ${u_{cf}}$ for ${k_\theta } \to \infty $ (i.e. clamped-free wing). This might appear as if the clamped-free wing is less stable than the pinned-free wing, which again looks counterintuitive. As shown in Fig. 2(d), for sufficiently large imperfections (i.e. very low values of ${k_\theta }$ ), the wing undergoes divergence at negligibly small flow velocities, and thus, flutter, if happens, would be the second instability; see the dashed line in Fig. 3, which shows the divergence speed boundary. On the other hand, for ${k_\theta } \to \infty $ , flutter may be the first or second instability, depending on the value of $\mu $ . As seen from the curves, ${u_{cf}}$ may undergo some sharp changes in the range of ${k_\theta } \in \left( {{{10}^{ - 2}},{{10}^2}} \right)$ . This is particularly visible for the lower mass ratios. These abrupt changes may be attributed to the ‘mode exchange’ or ‘role reversal’ phenomenon [Reference Païdoussis52], in which the unstable mode shifts from a lower mode to a higher mode or vice versa as a system parameter, such as ${k_\theta }$ , is varied. This can result in sharp variations in the critical flow velocity, even if the system parameter changes only slightly. For example, considering Fig. 4(a, b), the mode undergoing flutter shifts from second torsional mode ( ${u_{cf}} \simeq 2.02)$ to first torsional mode ( ${u_{cf}} \simeq 0.57$ ), while ${k_\theta }$ is increased from 1.2 to 1.3. Interestingly, within a narrow range of ${k_\theta }$ inside $\left( {{{10}^{ - 2}},{{10}^2}} \right)$ , by increasing ${k_\theta }$ , ${u_{cf}}$ decreases – the wing becomes dynamically less stable by decreasing support imperfections. A similar behaviour has previously been observed for imperfectly supported pipes conveying fluid [Reference Kheiri, Païdoussis, Costa Del Pozo and Amabili41].

Figure 4. Argand diagrams showing the mode-exchange or role reversal phenomenon: (a) ${{\rm{k}}_{\rm{\theta }}} = 1.2$ (second torsional mode flutter), and (b) ${{\rm{k}}_{\rm{\theta }}} = 1.3$ (first torsional mode flutter). The system parameters are the same as those in Table 1.

Figure 5 shows the variation of both ${u_{cf}}$ and ${u_{cd}}$ as a function of ${k_\theta }$ for several values of stiffness ratio. As seen from Fig. 5(a), for ${k_\theta } \gt 4.4$ , by increasing ${\rm{\Gamma }}$ , ${u_{cf}}$ decreases. Considering a fixed value for $EI$ , this means a wing with a lower $GJ$ (thus higher ${\rm{\Gamma }}$ ) becomes unstable at a lower flow velocity. This seems reasonable since for the configuration considered in this study, flutter occurs predominantly by a torsional mode. However, interestingly, a similar trend is not observed for lower values of ${k_\theta }$ . This is primarily due to the occurrence of the role reversal phenomenon. As seen from Fig. 5(b), for sufficiently low values of ${k_\theta }$ , regardless of the value of ${\rm{\Gamma }}$ , ${u_{cd}}$ approaches zero. On the other hand, for higher values of ${k_\theta }$ , ${u_{cd}}$ becomes non-zero, and it decreases as ${\rm{\Gamma }}$ is increased. This indicates that the divergence of the imperfectly supported wing is also predominantly influenced by the wing torsional stiffness. For a given ${k_\theta }$ and ${\rm{\Gamma }}$ , by comparing Fig. 5(a) and 5(b), one can find out whether flutter or divergence occurs first.

Figure 5. Variation of the dimensionless critical flow velocity for (a) flutter ( ${{\rm{u}}_{{\rm{cf}}}}$ ), and (b) divergence ( ${{\rm{u}}_{{\rm{cd}}}}$ ) of an imperfectly supported wing as a function of the end-spring stiffness ( ${{\rm{k}}_{\rm{\theta }}}$ ) for different values of stiffness ratio. The rest of the system parameters are the same as the dimensionless variables given in Table 1.

To further clarify the trend observed from Fig. 5(a) in the range ${k_\theta } \lt 4.4$ , Fig. 6 shows the Argand diagrams for a wing with ${k_\theta } = 1$ and for various values of ${\rm{\Gamma }}$ . As seen from the diagrams, flutter can occur not only via different torsional modes, it may also occur via bending modes. This observation further supports the frequent occurrence of the mode exchange (or role reversal) phenomenon in the present dynamical system. Another factor that should be taken into account is that except for the case with ${\rm{\Gamma }} = 1$ (Fig. 6(a)), divergence occurs first, and flutter occurs next, at higher flow velocities.

Figure 6. Argand diagrams showing the mode-exchange or role reversal phenomenon for ${{\rm{k}}_{\rm{\theta }}} = 1$ : (a) ${\rm{\Gamma }} = 1$ (first torsional mode flutter), (b) ${\rm{\Gamma }} = 2$ (second torsional mode flutter), (c) ${\rm{\Gamma }} = 3$ (second bending mode flutter), and (d) ${\rm{\Gamma }} = 4$ (first torsional mode flutter). The rest of the system parameters are the same as those in Table 1.

From Fig. 7, it is seen that, depending on the degree of imperfections, increasing ${r^2}$ may decrease or increase ${u_{cf}}$ . Similar mechanisms to those discussed in relation to Figs. 3 and 5 also operate here. Since ${r^2}$ only appears in the dynamics, it does not influence ${u_{cd}}$ , which is confirmed by a single dashed line shown in the figure.

Figure 7. Variation of the dimensionless critical flow velocity for flutter ( ${{\rm{u}}_{{\rm{cf}}}}$ ) of an imperfectly supported wing as a function of the end-spring stiffness ( ${{\rm{k}}_{\rm{\theta }}}$ ) for different values of dimensionless radius of gyration. Also, the dashed line shows the critical flow velocity for divergence, which is independent of ${{\rm{r}}^2}$ . The rest of the system parameters are the same as the dimensionless variables given in Table 1.

5.0 Concluding Remarks

The present paper, to the best of the authors’ knowledge, provided the first theoretical examination of the effects of support imperfections on the aeroelastic stability of uniform flexible wings. Following the extended Hamilton’s principle, the linear equations of motion were obtained in the dimensionless form, discretised spatially using Galerkin’s method, then cast in the state-space form and finally solved numerically to obtain the eigenvalues.

The numerical results showed that support imperfections can generally change the onset of the aeroelastic instabilities, their sequence, and the unstable mode, when compared to the wing with no support imperfections. The significance of the changes was found to be dependent on how large the imperfections were and to a lesser extent on the system parameters, such as mass ratio and stiffness ratio. For small support imperfections, small or no changes in the aeroelastic behaviour were observed. However, for large support imperfections, i.e. loosely supported wings, dramatic and sometimes unexpected changes may take place. For example, it was found that because of the role reversal or mode exchange phenomenon, flutter speed may change sharply when the support stiffness is varied slightly within a narrow range. Also, in some scenarios, by decreasing imperfections (i.e., stiffer support), the critical flow velocity for flutter is decreased.

Nevertheless, the linear theory predictions cannot be trusted beyond the first instability, and thus, a nonlinear model should be sought in future to verify the linear model. In addition, experimental studies should be conducted to validate the numerical results presented in this paper. The authors hope that this first study can stimulate further research on the topic in the community.

Acknowledgements

The support for this research project provided by the Natural Sciences and Engineering Research Council of Canada (NSERC), as well as Concordia University, is acknowledged.

Competing interests

The authors declare no competing interests.

Appendix

6.0 First-order matrices

The matrices ${\mathcal A}$ , ${\mathcal B}$ and vector ${\mathcal F}$ from the first-order equations of motion are given as follows:

(19) \begin{align}{\mathcal A} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}{{{\mathbf{M}}^{\left( 1 \right)}}}& {}{{{\mathbf{M}}^{\left( 2 \right)}}} {}&{{{\mathbf{C}}^{\left( 1 \right)}}}& {}{{{\mathbf{C}}^{\left( 2 \right)}}} {}&0 {}&0 {}&0 {}&0\\{{{\mathbf{M}}^{\left( 3 \right)}}}& {}{{{\mathbf{M}}^{\left( 4 \right)}}} {}&{{{\mathbf{C}}^{\left( 3 \right)}}} {}&{{{\mathbf{C}}^{\left( 4 \right)}}} {}&0 {}&0 {}&0& {}0\\0 {}&0 {}&{\mathbf{I}} {}&0 {}&0 {}&0 {}&0 {}&0\\0 {}&0 {}&0 {}&{\mathbf{I}} {}&0 {}&0 {}&0 {}&0\\0 {}&0 {}&0 {}&0 {}&{\mathbf{I}} {}&0 {}&0 {}&0\\0 {}&0 {}&0 {}&0 {}&0 {}&{\mathbf{I}} {}&0 {}&0\\0 {}&0 {}&0 {}&0 {}&0 {}&0 {}&{\mathbf{I}} {}&0\\0 {}&0 {}&0 {}&0 {}&0 {}&0 {}&0 {}&{\mathbf{I}}\end{array}} \right]\end{align}
(20) \begin{align}{\mathcal B} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}0 {}&0 {}&{ - {{\mathbf{K}}^{\left( 1 \right)}}} {}&{ - {{\mathbf{K}}^{\left( 2 \right)}}} {}&{ - {{\mathbf{G}}^{\left( 1 \right)}}} {}&{ - {{\mathbf{G}}^{\left( 2 \right)}}} {}&{ - {{\mathbf{G}}^{\left( 3 \right)}}} {}&{ - {{\mathbf{G}}^{\left( 4 \right)}}}\\0& {}0& {}{ - {{\mathbf{K}}^{\left( 3 \right)}}}& {}{ - {{\mathbf{K}}^{\left( 4 \right)}}}& {}{ - {{\mathbf{G}}^{\left( 5 \right)}}}& {}{ - {{\mathbf{G}}^{\left( 6 \right)}}}& {}{ - {{\mathbf{G}}^{\left( 7 \right)}}}& {}{ - {{\mathbf{G}}^{\left( 8 \right)}}}\\{\mathbf{I}}& {}0& {}0 {}&0 {}&0& {}0 {}&0& {}0\\0 {}&{\mathbf{I}} {}&0 {}&0 {}&0 {}&0 {}&0 {}&0\\0 {}&0& {}{\mathbf{I}}& {}0& {}{ - {\varepsilon _1}{\mathbf{I}}}& {}0& {}0& {}0\\0 {}&0& {}{\mathbf{I}}& {}0& {}0& {}{ - {\varepsilon _2}{\mathbf{I}}}& {}0& {}0\\0 {}&0& {}0& {}{\mathbf{I}}& {}0 {}&0& {}{ - {\varepsilon _1}{\mathbf{I}}}& {}0\\0 {}&0 {}&0& {}{\mathbf{I}}& {}0 {}&0& {}0& {}{ - {\varepsilon _2}{\mathbf{I}}}\end{array}} \right]\end{align}
(21) \begin{align}{\mathcal F} = \left[ {\begin{array}{*{20}{c}}{ - {{\mathbf{F}}^{\left( 1 \right)}}}\\{ - {{\mathbf{F}}^{\left( 2 \right)}}}\\0\\0\\0\\0\\0\\0\end{array}} \right]{\rm{\;\;\;\;}},\end{align}

where $0$ and ${\bf{I}}$ are $N$ -by- $N$ zero and identity matrices, respectively.

Footnotes

1 Similarly, ${({\rm{\;\;}})^{{\rm{\prime\prime}}}} = {\partial ^2}\left( {{\rm{\;\;}}} \right)/\partial {y^2}$ and ${({\rm{\;\;}})^{{\rm{\prime\prime\prime\prime}}}} = {\partial ^4}\left( {{\rm{\;\;}}} \right)/\partial {y^4}$ .

References

Bisplinghoff, R.L., Ashley, H. and Halfman, R.L. Principles of Aeroelasticity, New York: Dover Publications, Inc., 1996.Google Scholar
Hodges, D.H. and Pierce, G.A. Introduction to Structural Dynamics and Aeroelasticity, 2nd ed, New York, NY, USA: Cambridge University Press, 2011.10.1017/CBO9780511997112CrossRefGoogle Scholar
Federal Aviation Administration, Part 25—Airworthiness Standards: Transport Category Airplanes, USA, 2025. https://www.ecfr.gov/current/title-14/part-25. Accessed 16, April, 2025.Google Scholar
Peters, D.A., Karunamoorthy, S. and Caom, W.-M. Finite state induced flow models. i-two-dimensional thin airfoil, J. Aircr., 1995, 32, (2), pp 313322.10.2514/3.46718CrossRefGoogle Scholar
Afonso, F., Vale, J., Oliveira, É, Lau, F. and Suleman, A. A review on non-linear aeroelasticity of high aspect-ratio wings, Prog. Aerosp. Sci., 2017, 89, pp 4057.10.1016/j.paerosci.2016.12.004CrossRefGoogle Scholar
Tang, D. and Dowell, E.H. Experimental and theoretical study on aeroelastic response of high-aspect-ratio wings, AIAA J., 2001, 39, (8), pp 14301441.10.2514/2.1484CrossRefGoogle Scholar
Patil, M.J., Hodges, D.H. and Cesnik, C.E.S. Nonlinear aeroelasticity and flight dynamics of high-altitude long-endurance aircraft, J. Aircr., 2001, 38, (1), pp 8894.10.2514/2.2738CrossRefGoogle Scholar
Qin, Z. and Librescu, L. Dynamic aeroelastic response of aircraft wings modeled as anisotropic thin-walled beams, J. Aircr., 2003, 40, (3), pp 532543.10.2514/2.3127CrossRefGoogle Scholar
Tang, D.M. and Dowell, E.H. Effects of geometric structural nonlinearity on flutter and limit cycle oscillations of high-aspect-ratio wings, J. Fluids Struct., 2004, 19, (3), pp 291306.10.1016/j.jfluidstructs.2003.10.007CrossRefGoogle Scholar
Shams, S., Lahidjani, M.S. and Haddadpour, H. Nonlinear aeroelastic response of slender wings based on Wagner function, Thin-Walled Struct., 2008, 46, (11), pp 11921203.10.1016/j.tws.2008.03.001CrossRefGoogle Scholar
Berci, M. and Cavallaro, R. A hybrid reduced-order model for the aeroelastic analysis of flexible subsonic wings—a parametric assessment, Aerospace, 2018, 5, (3), p 76.10.3390/aerospace5030076CrossRefGoogle Scholar
Amoozgar, M., Fazelzadeh, S., Khodaparast, H.H., Friswell, M. and Cooper, J. Aeroelastic stability analysis of aircraft wings with initial curvature, Aerosp. Sci. Technol., 2020, 107, p 106241.10.1016/j.ast.2020.106241CrossRefGoogle Scholar
Zhang, X., Kheiri, M. and Xie, W.-F. Nonlinear dynamics and gust response of a two-dimensional wing, Int. J. Non-Linear Mech., 2020, 123, p 103478.10.1016/j.ijnonlinmec.2020.103478CrossRefGoogle Scholar
Zhang, X., Kheiri, M. and Xie, W.-F. Aeroservoelasticity of an airfoil with parametric uncertainty and subjected to atmospheric gusts, AIAA J., 2021, 59, (11), pp 43264341.10.2514/1.J060089CrossRefGoogle Scholar
Lee, B., Price, S. and Wong, Y. Nonlinear aeroelastic analysis of airfoils: bifurcation and chaos, Prog. Aerosp. Sci., 1999, 35, (3), pp 205334.10.1016/S0376-0421(98)00015-3CrossRefGoogle Scholar
Li, D., Guo, S. and Xiang, J. Aeroelastic dynamic response and control of an airfoil section with control surface nonlinearities, J. Sound Vib., 2010, 329, (22), pp 47564771.10.1016/j.jsv.2010.06.006CrossRefGoogle Scholar
Sales, T.P., Marques, F.D., Pereira, D.A. and Rade, D.A. Dynamic assessment of nonlinear typical section aeroviscoelastic systems using fractional derivative-based viscoelastic model, J. Sound Vib., 2018, 423, pp 230245.10.1016/j.jsv.2018.02.008CrossRefGoogle Scholar
Panchal, J. and Benaroya, H. Review of control surface freeplay, Prog. Aerosp. Sci., 2021, 127, p 100729.10.1016/j.paerosci.2021.100729CrossRefGoogle Scholar
Bueno, D.D., Wayhs-Lopes, L.D. and Dowell, E.H. Control-surface structural nonlinearities in aeroelasticity: a state of the art review, AIAA J., 2022, 60 (6), pp 33643376.10.2514/1.J060714CrossRefGoogle Scholar
da Silva, J.A.I. and Marques, F.D. Characterization of typical aeroelastic sections under combined structural concentrated nonlinearities, J. Vib. Control, 2022, 28, (13–14), pp 17391753.10.1177/10775463211000161CrossRefGoogle Scholar
Amoozgar, M., Castrichini, A., Garvey, S., Friswell, M., Cooper, J. and Ajaj, R. The effect of a nonlinear energy sink on the gust response of a wing, Aerosp. Sci. Technol., 2024, 145, p 108904.10.1016/j.ast.2024.108904CrossRefGoogle Scholar
Yuan, W., Sandhu, R. and Poirel, D. Fully coupled aeroelastic analyses of wing flutter towards application to complex aircraft configurations, J. Aerosp. Eng., 2021, 34, (2), p 04020117.10.1061/(ASCE)AS.1943-5525.0001232CrossRefGoogle Scholar
Farhat, C. and Lesoinne, M. Two efficient staggered algorithms for the serial and parallel solution of three-dimensional nonlinear transient aeroelastic problems, Comput. Methods Appl. Mech. Eng., 2000, 182, (3–4), pp 499515.10.1016/S0045-7825(99)00206-6CrossRefGoogle Scholar
Farhat, C., Geuzaine, P. and Brown, G. Application of a three-field nonlinear fluid–structure formulation to the prediction of the aeroelastic parameters of an f-16 fighter, Comput. Fluids, 2003, 32, (1), pp 329.10.1016/S0045-7930(01)00104-9CrossRefGoogle Scholar
Dang, H., Yang, Z. and Li, Y. Accelerated loosely-coupled cfd/csd method for nonlinear static aeroelasticity analysis, Aerosp. Sci. Technol., 2010, 14, (4), pp 250258.10.1016/j.ast.2010.01.004CrossRefGoogle Scholar
Albano, E. and Rodden, W.P. A doublet-lattice method for calculating lift distributions on oscillating surfaces in subsonic flows, AIAA J., 1969, 7 (2), pp 279285.10.2514/3.5086CrossRefGoogle Scholar
Katz, J. and Plotkin, A. Low-Speed Aerodynamics, Vol. 13, New York, NY, USA: Cambridge University Press, 2001.10.1017/CBO9780511810329CrossRefGoogle Scholar
Murua, J., Palacios, R. and Graham, J.M.R. Applications of the unsteady vortex-lattice method in aircraft aeroelasticity and flight dynamics, Prog. Aerosp. Sci., 2012, 55, pp 4672.10.1016/j.paerosci.2012.06.001CrossRefGoogle Scholar
ZONA Technology, Inc. Tech. Rep. ZAERO Version 9.3, User’s manual, Scottsdale, AZ, USA, 2020. https://www.zonatech.com/zaero.html/. Accessed 15, April, 2025.Google Scholar
Parenteau, M. and Laurendeau, E. A general modal frequency-domain vortex lattice method for aeroelastic analyses, J. Fluids Struct., 2020, 99, p 103146.10.1016/j.jfluidstructs.2020.103146CrossRefGoogle Scholar
Hilger, J. and Ritter, M.R. Nonlinear aeroelastic simulations and stability analysis of the pazy wing aeroelastic benchmark, Aerospace, 2021, 8, (10), p 308.10.3390/aerospace8100308CrossRefGoogle Scholar
Meher-Homji, C.B., Gabriles, G., et al., Gas turbine blade failures-causes, avoidance, and troubleshooting, in The 27th Turbomachinery Symposium, Texas, USA: Texas A&M University, Turbomachinery Laboratories, 1998, pp 129180.Google Scholar
Munns, T.E. and Kent, R.M. Structural health monitoring: degradation mechanisms and system requirements, in 19th Digital Avionics Systems Conference. Proceedings (Cat. No. 00CH37126), Vol. 2, IEEE, 2000, Philadelphia, PA, USA, pp 6C2-18.10.1109/DASC.2000.884907CrossRefGoogle Scholar
Wang, K. and Inman, D.J. Crack-induced effects on aeroelasticity of an unswept composite wing, AIAA J., 2007, 45, (3), pp 542551.10.2514/1.21689CrossRefGoogle Scholar
Hoseini, H.S. and Hodges, D.H. Aeroelastic stability analysis of damaged high-aspect-ratio composite wings, J. Aircr., 2019, 56, (5), pp 17941808.10.2514/1.C035098CrossRefGoogle Scholar
Torabi, A., Shams, S., Narab, M.F. and Atashgah, M.A. Unsteady aero-elastic analysis of a composite wing containing an edge crack, Aerosp. Sci. Technol., 2021, 115, p 106769.10.1016/j.ast.2021.106769CrossRefGoogle Scholar
Torabi, A.R., Shams, S., Fatehi Narab, M. and Atashgah, M.A. Delamination effects on the unsteady aero-elastic behavior of composite wing by modal analysis, J. Vib. Control, 2022, 28 (19–20), pp 27322745.10.1177/10775463211019213CrossRefGoogle Scholar
Kheiri, M. and Païdoussis, M.P. Dynamics and stability of a flexible pinned-free cylinder in axial flow, J. Fluids Struct., 2015, 55, pp 204217.10.1016/j.jfluidstructs.2015.02.013CrossRefGoogle Scholar
Tabatabaei, S.S., Kheiri, M. and Dargahi, J. Dynamics and stability of imperfect flexible cylinders in axial flow, J. Fluids Struct., 2021, 105, p 103321.10.1016/j.jfluidstructs.2021.103321CrossRefGoogle Scholar
Kheiri, M., Païdoussis, M.P., Costa Del Pozo, G. and Amabili, M. Dynamics of a pipe conveying fluid flexibly restrained at the ends, J. Fluids Struct., 2014, 49, pp 360385.10.1016/j.jfluidstructs.2013.11.023CrossRefGoogle Scholar
Kheiri, M. Nonlinear dynamics of imperfectly-supported pipes conveying fluid, J. Fluids Struct., 2020, 93, p 102850.10.1016/j.jfluidstructs.2019.102850CrossRefGoogle Scholar
Riazat, M. and Kheiri, M. Three-dimensional nonlinear dynamics of imperfectly supported pipes conveying fluid, J. Fluids Struct., 2023, 123, p 104011.10.1016/j.jfluidstructs.2023.104011CrossRefGoogle Scholar
Asadi, D. and Farsadi, T. Active flutter control of thin walled wing-engine system using piezoelectric actuators, Aerosp. Sci. Technol., 2020, 102, p 105853.10.1016/j.ast.2020.105853CrossRefGoogle Scholar
Wang, X., Zhou, W., Mu, R. and Wu, Z. A new deformation control approach for flexible wings using moving masses, Aerosp. Sci. Technol., 2020, 106, p 106118.10.1016/j.ast.2020.106118CrossRefGoogle Scholar
Silva, G.C., Donadon, M.V. and Silvestre, F.J. Experimental and numerical investigations on the nonlinear aeroelastic behavior of high aspect-ratio wings for different chord-wise store positions under stall and follower aerodynamic load models, Int. J. Non-Linear Mech., 2021, 131, p 103685.10.1016/j.ijnonlinmec.2021.103685CrossRefGoogle Scholar
Vindigni, C.R., Mantegna, G., Orlando, C. and Alaimo, A. Simple adaptive wing-aileron flutter suppression system, J. Sound Vib., 2024, 570, p 118151.10.1016/j.jsv.2023.118151CrossRefGoogle Scholar
Modaress-Aval, A.H., Bakhtiari-Nejad, F., Dowell, E.H., Peters, D.A. and Shahverdi, H. A comparative study of nonlinear aeroelastic models for high aspect ratio wings, J. Fluids Struct., 2019, 85, pp 249274.10.1016/j.jfluidstructs.2019.01.003CrossRefGoogle Scholar
Computational Fluid Dynamics Committee, Guide for the verification and validation of computational fluid dynamics simulations (AIAA G-077-1998) (1998).Google Scholar
Goland, M. The flutter of a uniform cantilever wing, J. Appl. Mech., 1945, 12 (4), pp A197A208.10.1115/1.4009489CrossRefGoogle Scholar
Haddadpour, H. and Firouz-Abadi, R. Evaluation of quasi-steady aerodynamic modeling for flutter prediction of aircraft wings in incompressible flow, Thin-Walled Struct., 2006, 44 (9), pp 931936.10.1016/j.tws.2006.08.020CrossRefGoogle Scholar
Païdoussis, M.P. Fluid-Structure Interactions: Slender Structures and Axial Flow. Volume 1, 2nd ed, London: Academic Press, 2014.Google Scholar
Figure 0

Figure 1. Schematic drawing of a uniform flexible wing imperfectly supported at ${\rm{y}} = 0$ and free at ${\rm{y}} = \ell $, where $\ell $ is the span. Two rotational springs of the stiffness ${{\rm{K}}_{\rm{w}}}$ and ${{\rm{K}}_{\rm{\theta }}}$ are used to model the support imperfection; the springs are attached to the wing at one end and attached to a rigid support (e.g., fuselage) at the other end; ${\rm{xyz}}$ is the Cartesian coordinate system attached to the undeformed wing; also, ${\rm{c}}$ is the chord length and ${\rm{U}}$ is the freestream velocity.

Figure 1

Table 1. Parameters of the wing used in Ref. [7] and the corresponding dimensionless parameters

Figure 2

Table 2. Numerical convergence study using a different number of mode shapes for a wing with the parameters given in Table 1; also, ${{\rm{k}}_{\rm{\theta }}} = {10^6}$

Figure 3

Table 3. Comparison between the critical speeds for the aeroelastic instabilities of two wings obtained from the present model and those found from the literature

Figure 4

Figure 2. Argand diagrams showing the evolution of the first few dynamical modes of an imperfectly-supported wing with parameters given in Table 1 for: (a) ${{\rm{k}}_{\rm{\theta }}} = {10^6}$, (b) ${{\rm{k}}_{\rm{\theta }}} = {10^0}$, (c) ${{\rm{k}}_{\rm{\theta }}} = {10^{ - 1}}$, and (d) ${{\rm{k}}_{\rm{\theta }}} = {10^{ - 6}}$. The numerals in the plots, close to the loci correspond to the values of the dimensionless flow velocity. The ${\rm{x}}$ and ${\rm{y}}$ axes have been normalised with respect to the natural frequency of the first torsional mode of a cantilevered beam; also, ${{\rm{u}}_{{\rm{cd}}}}$ and ${{\rm{u}}_{{\rm{cf}}}}$ represent the critical value of ${\rm{u}}$ for divergence and flutter, respectively.

Figure 5

Figure 3. Variation of the dimensionless critical flow velocity for flutter (${{\rm{u}}_{{\rm{cf}}}}$) of an imperfectly supported wing as a function of the end-spring stiffness (${{\rm{k}}_{\rm{\theta }}}$) for different values of mass ratio. Also, the dashed line shows the critical flow velocity for divergence, which is independent of ${\rm{\mu }}$. The rest of system parameters are the same as the dimensionless variables given in Table 1.

Figure 6

Figure 4. Argand diagrams showing the mode-exchange or role reversal phenomenon: (a) ${{\rm{k}}_{\rm{\theta }}} = 1.2$ (second torsional mode flutter), and (b) ${{\rm{k}}_{\rm{\theta }}} = 1.3$ (first torsional mode flutter). The system parameters are the same as those in Table 1.

Figure 7

Figure 5. Variation of the dimensionless critical flow velocity for (a) flutter (${{\rm{u}}_{{\rm{cf}}}}$), and (b) divergence (${{\rm{u}}_{{\rm{cd}}}}$) of an imperfectly supported wing as a function of the end-spring stiffness (${{\rm{k}}_{\rm{\theta }}}$) for different values of stiffness ratio. The rest of the system parameters are the same as the dimensionless variables given in Table 1.

Figure 8

Figure 6. Argand diagrams showing the mode-exchange or role reversal phenomenon for ${{\rm{k}}_{\rm{\theta }}} = 1$: (a) ${\rm{\Gamma }} = 1$ (first torsional mode flutter), (b) ${\rm{\Gamma }} = 2$ (second torsional mode flutter), (c) ${\rm{\Gamma }} = 3$ (second bending mode flutter), and (d) ${\rm{\Gamma }} = 4$ (first torsional mode flutter). The rest of the system parameters are the same as those in Table 1.

Figure 9

Figure 7. Variation of the dimensionless critical flow velocity for flutter (${{\rm{u}}_{{\rm{cf}}}}$) of an imperfectly supported wing as a function of the end-spring stiffness (${{\rm{k}}_{\rm{\theta }}}$) for different values of dimensionless radius of gyration. Also, the dashed line shows the critical flow velocity for divergence, which is independent of ${{\rm{r}}^2}$. The rest of the system parameters are the same as the dimensionless variables given in Table 1.