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 ¹

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 Farsadi44–Reference 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.