NOMENCLATURE
- AoA
-
Angle-of-Attack
- CAD
-
Computer-Aided Design
- CFD
-
Computational Fluid Dynamics
- DLM
-
Doublet Lattice Method
- FEM
-
Finite Element Model
- GVT
-
Ground Vibration Test
- HALE
-
High-Altitude Long-Endurance
- HARW
-
High-Aspect-Ratio Wing
- LCO
-
Limit Cycle Oscillation
- L/D
-
Lift-to-Drag Ratio
- MALE
-
Medium-Altitude Long-Endurance
- MIL STD
-
Military Standard
- MTBL
-
Mean Time Before Loss
- MTBCF
-
Mean Time before Critical Failure
- ROM
-
Reduced-Order Model
- UAV
-
Unmanned Aerial Vehicle
Latin symbols
- b
-
semi-span of wing
-
$C_{L\alpha}$
-
lift curve slope
- D
-
damping terms
- e
-
eccentricity
- h
-
plunge motion
-
$I_{\alpha}$
-
moment of inertia
- k
-
reduced frequency
-
$k_h$
-
bending stiffness
-
$k_\alpha$
-
torsional stiffness
- L
-
aerodynamic lift
- m
-
mass
- M
-
aerodynamic moment
- q
-
dynamic pressure
-
$q_{12}$
-
generalised coordinates
- Q
-
generalised force
-
$S_{\alpha}$
-
mass unbalance
- T
-
kinetic energy
- u
-
displacement in transverse direction
- U
-
potential energy
- V
-
velocity
- w
-
rotation about longitudinal axis
- W
-
work done
Greek symbols
-
$\alpha$
-
angle-of-attack (pitching motion)
-
$\Gamma$
-
real part of eigenvalue
-
$\gamma$
-
rate of decay
-
$\lambda$
-
eigenvalues
-
$\sigma$
-
shifting function
-
$\Omega$
-
imaginary part of eigenvalue
1.0 INTRODUCTION
UAVs are becoming increasingly common due to their ease of deployment, low cost of maintenance, manufacturing, operation, platform and ground control and good mobility. UAVs offer longer endurance and high reliability by having shorter mean time before loss (MTBL), and mean time before critical failure (MTBCF)(Reference Tsach, Tatievsky and London1). UAVs are commonly used for operations that are regarded as dull, dirty and dangerous. This refers to long and exhausting missions posing a high risk to operators. Acceptance of the use of UAVs grew when the benefits of this platform were realised, being complementary to piloted missions. Apart from military applications, civil infrastructure alone is expected to dominate the market of more than
$\$$
45 Bn for UAV usage in the near future(Reference Shakhatreh, Sawalmeh, Al-Fuqaha, Dou, ALmaita, Khalil, Othman, Shamsiah, Khreishah and Guizani2). MALE and HALE vehicles in particular are evolving as two main categories of unmanned aerial platform. UAVs are extensively used for city and coastal surveillance, precision agriculture, infrastructure inspection, wireless coverage, goods delivery and traffic monitoring(Reference Mohammad, Idries, Mohamed, Al-JAroodi and Jwahar3,Reference Gomez-Candon, De Castro and Lopez-Granados4) . Military applications of these UAVs include surveillance and treaty monitoring, assessment of hostile areas and search and rescue (SAR) operations(Reference Panagiotou, Tsavlidis and Yakinth5).
The endurance of UAVs can be increased by tailoring the aerofoil and wing geometries to achieve lower drag in dominant operating regimes. High-endurance UAVs are characterised by long, slender wings, as most UAV platforms are designed to operate in slow, subsonic regimes that offer performance benefits for a high-aspect-ratio wing (HARW)(Reference Gorniak, Goraj and Olszanski6). The lift-to-drag ratio (L/D) is a key performance parameter that directly affects the efficiency of an aircraft(Reference Anderson7,Reference Raymer8,Reference Panagiotou, Kaparos and Yakinthos9) . Higher efficiency means lower fuel consumption and hence less environmental pollution and greater endurance(Reference Tsach, Tatievsky and London1). HARWs offer an improved L/D ratio by reducing the induced drag, which accounts for almost half of the total drag experienced by aircraft flying under cruise conditions(Reference Abbas, De Vicente and Valero10). However, the use of HARWs is associated with a trade-off with structural efficiency and safety. Longer wings tend to be more flexible and deform easily under load, hence being more susceptible to detrimental aeroelastic effects(Reference Afonso, Vale, Oliveira, Lau and Suleman11). Divergence, control reversal and flutter are some major aeroelastic effects that can range from mere discomfort to complete destruction of the body in flight (Fig. 1)(Reference Kehoe12).
It is these aeroelastic phenomena that mostly constrain the flight envelope, or more importantly the maximum flight speed, of any vehicle (Fig. 2)(Reference Anderson7). Aeroelastic effects can result in divergence, control reversal and flutter. Theoretically, divergence occurs at the point where the aeroelastic twist of a structure becomes infinite and the structure fails(Reference Dowell14). In control reversal, the effectiveness of the control surfaces, such as ailerons, is completely reduced or even reversed from their normal operation(Reference Fung15). Flutter is a dynamic aeroelastic phenomenon in which the aerodynamic, elastic and inertial forces interact. In this effect, energy from the flow is rapidly absorbed and transformed into mechanical vibrations that result in a diverging vibrational response.
Figure 1. Venn diagram showing all possible relations between the set of disciplines that collectively form the basis of aeroelasticity (adopted from Hodges)(Reference Hodges and Pierce13).
Figure 2. V–n diagram for an aircraft, showing the load factor (G) boundaries that should not be exceeded during flight. The aerodynamic limit marks the stall region, whereas the load limit defines the structural stress limits. The maximum allowed airspeed to avoid aeroelastic effects is shown as the airspeed limit.
Aeroelasticity is an inherent property of the structure and cannot be eliminated(Reference Shubov16). However, it must be avoided either by the design of the structure or the use of some active/passive control mechanism(Reference Anderson, White and Neal17–Reference Ghasemikaram, Mazidi, Fazel and Fazelzadeh21) that prevents or at least reduces its harmful effects. While the effects of aeroelasticity are experienced over the complete body of an aircraft, its effects are strongest on the wings and tail. These areas are highly prone to failure, as these are subject to such effects early in the flight regime.
The concept of aeroelasticity is as old as the first aircraft, or even older than that(Reference Dowell14). Various methods such as mathematical modelling, GVTs, bending tests, wind-tunnel experiments and flight flutter testing constitute an essential phase of the aircraft design and development process(Reference Friedmann22,Reference Garrick and Reed23) . Most aeroelastic investigations have focused on flutter calculations, because a century of experience suggests that the onset of flutter for a conventional wing appears well before other aeroelastic effects(Reference Edwards and Wieseman24). Most modern techniques are based on a framework of divergence analysis within the confines of the flutter methodology. While methods for flutter suppression are abundantly available in literature(Reference Anderson, White and Neal17,Reference Livne25) , it is imperative to first establish the natural aeroelastic characteristics of the wing being designed.
A good design that is developed considering the intended operating conditions does not experience flutter within its flight envelope. Some researches(Reference Kehoe12) have aimed to devise improved models, approaches and method for such testing, as well as techniques for the suppression and control of aeroelastic effects, as reported abundantly in literature(Reference Njuguna26,Reference Schuster, Liu and Huttsell27) . Zhang et al.(Reference Zhang, Wang, Ye and Quan28) constructed a non-linear unsteady aerodynamic ROM using a Radial Basis Function (RBF) neural network model to investigate the LCOs. The time-marching solutions obtained through a hybrid linear multi-step algorithm agreed well with those obtained from CFD, with a reduction in computational time by one to two orders. Eller(Reference Eller29) presented a robust method by defining flutter as an eigenvalue problem of piecewise-quadratic nature, which yielded results within reasonable agreement with the pk-method. This approach overcame the inability of other methods to include damping using data corresponding to pure harmonic motion. Recently, Voss et al.(Reference Voss, Schaefer and Vidy30) carried out an extensive numerical study on the flutter characteristics of a DLR-F19/SACCON flying wing configuration. The results obtained from the classical DLM as well as CFD-based methods were compared, revealing close agreement regarding the flutter mechanisms. However, variations in the flutter speed were observed, depending upon the aerodynamic formulation used. Jacek et al.(Reference Gadomski, Hernik and Goraj31) proposed an effective design methodology for loaded structures of UAVs, especially MALE/HALE platforms, to obtain designs that are well suited to their intended mission capability. The insufficient torsional rigidity of their design was addressed by joining additional layers in a sandwich structure for the skin, thus demonstrating the ease of aeroelastic tailoring in composite structures.
We present herein a novel, semi-analytical approach to examine the flutter characteristics of the HARW of a MALE UAV. The construction of the wing demonstrates those used in modern flying platforms, integrating metals with composite structures. The complex task of flutter modelling is made less complicated by converting the detailed computer-aided design (CAD) model into ROMs, which is widely recognised to be a computationally efficient approach(Reference Amsallem and Farhat32). The flutter response is calculated using the analytical p-method, where the stiffness as well as geometric and modal parameters are incorporated from the FE model. Subsequently, the NX Nastran solver is used to computationally ascertain the stability response of the wing at a range of altitudes using the k- and pk-methods. An extended eigenvalue analysis is applied to determine the flutter dynamic pressures and velocities, which define the stability boundary of the wing at each altitude. The work is further extended to include a parametric study, which investigates the effect of the geometric properties of the wing on the flutter speed. This provides a viable approach that can be adopted for optimisation of flutter speeds by varying the thickness of structural members. The work begins with the derivation of the two-degree-of-freedom (2-DOF) pitch and plunge motion, which allows for bending and torsional deformation. The flutter framework derived from the Lagrange equations leads to the formulation of the flutter eigenvalue problem. The eigenvalue solution, though robust, yields only qualitative information about the stability of the system at defined flight conditions. However, it does not specify the stability boundaries, which are the design objective in this work. Therefore, such analysis is augmented by extended eigenvalue analysis, yielding the critical flutter velocities. The geometric and structural parametric information required for this analysis are computed computationally using the FE model. The frequencies and shapes for the first ten modes are calculated for use in the subsequent flutter analysis. The analytical results are validated numerically using simulations in Siemen’s FEMAP (NX Nastran) software with the V–g approach, revealing close agreement between the results of the p-, k-, and pk-methods. The effect of varying the thickness of the ribs and spars on the aeroelastic response is considered, bearing in mind the trade-off between the flutter speed and structural weight.
The remainder of this paper is organised as follows: Section 2 presents the analytical modelling of the wing flutter system to derive the governing equations. The finite element modelling and computational analysis approach are outlined in Section 3, followed by the modal analysis in Section 4. The results of the analytical and computational analyses are discussed and compared in Section 5. The procedure for and results of the parametric analysis are covered in Section 6. A commentary on the findings of the whole study, along with recommendations for future analysis, are given in Section 7.
2.0 ANALYTICAL MODELLING
Flutter occurs due to the coalescence of any two of the many natural modes of vibration. However, it is reasonable to consider only the first few modes when evaluating aeroelastic stability characteristics. Flutter can be analysed by treating the wing as a 2-DOF system exhibiting plunge and pitch motions, which correspond to span-wise bending and torsion. To define the 2-DOF plunge and pitch flutter model, consider the typical section of a wing shown in Fig. 3. To derive the equations of motion, Lagrange’s principle is used based on the total Kinetic Energy (KE) and Potential Energy (PE) of the system. The first step of this process requires a representation of the pitch and plunge motion in generalised coordinates (Reference Dowell14):
\begin{eqnarray}u = x[{\rm cos}\alpha-1] \cong 0, \end{eqnarray}
\begin{eqnarray}w = -h-x{\rm sin}\alpha \cong 0 \end{eqnarray}
Figure 3. A two-dimensional wing model with pitch (
$\alpha$
) and plunge (h) DOF. The elastic axis is offset from the aerodynamic centre and the centre of gravity by a distance e and x, respectively. (Adopted from Dowell.)(Reference Gomez-Candon, De Castro and Lopez-Granados4).
Once the displacement functions have been defined, they can be used to evaluate the expressions for the KE and PE. The bending and torsional stiffness are denoted as
$k_h$
and
$k_\alpha$
, respectively. The resulting expressions for the KE and PE can then be applied in Lagrange’s principle.
\begin{equation}T=\frac{1}{2} \dot{h}^{2} m+\frac{1}{2} 2 \dot{h} \dot{\alpha} S_{\alpha}+\frac{1}{2}\dot{\alpha}^{2} l_{\alpha}, \end{equation}
\begin{equation}U=\frac{1}{2} K_{h} h^{2}+\frac{1}{2} K_{\alpha} \alpha^{2} \end{equation}
Using the expression for the KE and PE in Lagrange’s principle:
\begin{equation} -\frac{{\rm d}}{{\rm d t}}\left(\frac{\partial(T-U)}{\partial \dot{h}}\right)+\frac{\partial(T-U)}{\partial h}+Q_{h}=0 \end{equation}
\begin{equation}-\frac{{\rm d}}{{\rm d t}}\left(\frac{\partial(T-U)}{\partial \dot{\alpha}}\right)+\frac{\partial(T-U)}{\partial \alpha}+Q_{\alpha}=0, \end{equation}
where
$Q_h$
and
$Q_\alpha$
are the generalised forces derived by evaluating the work done by aerodynamic forces, which is given by
\begin{eqnarray}\delta W=Q_{h} \delta h+Q_{\alpha} \delta \alpha \end{eqnarray}
The aerodynamics forces act as a resultant of the pressure distribution over the surface. This suggests that the lift and moment created by the aerodynamic forces are actually the generalised forces that correspond to the pitch and plunge motion(Reference Dowell14).
\begin{align}m \ddot{h}+K_{h} h+S_{\alpha} \ddot{\alpha}=-L,\end{align}
\begin{align}S_{\alpha} \ddot{h}+I_{\alpha} \ddot{\alpha}+K_{\alpha} \alpha={\rm M}_{y}\end{align}
Many aerodynamic theories available in literature can be used in Equations (8) and (9) to describe the aerodynamic lift and moment. For the geometric configuration of a wing alone, classical, steady aerodynamic theory can be applied.
\begin{equation}L=q S C_{L \alpha} \alpha, \end{equation}
\begin{equation}{\rm M}_{y}=e L \end{equation}
In this way, Equations (8) and (9) can be written in the form
\begin{equation}m \ddot{h}+S_{\alpha} \ddot{\alpha}+K_{h} h+q S C_{L \alpha} \alpha=0, \end{equation}
\begin{equation}I_{\alpha} \ddot{\alpha}+S_{\alpha} \ddot{h}+K_{\alpha} \alpha-q S e C_{L \alpha}=0 \end{equation}
Simple harmonic motion is assumed for the stability analysis, stated mathematically as
\begin{equation}h=\bar{h} e^{p t}, \end{equation}
\begin{equation}\alpha=\bar{\alpha} e^{p t}, \end{equation}
Here, p is a complex number whose real part determines the converging or diverging nature of the motion. Substitution of h and
$\alpha$
into Equations (12) and (13) results in a set of equations that can be written in matrix form as
\begin{eqnarray}\left[\begin{array}{cc}m p^{2}+K_{h} & S_{\alpha} p^{2}+q S C_{L \alpha} \\[3pt] S_{\alpha} p^{2} & I_{\alpha} p^{2}+k_{\alpha}-q S e C_{L \alpha}\end{array}\right]\left\{\begin{array}{l}\bar{h} e^{p t} \\[3pt] \bar{\alpha} e^{p t}\end{array}\right\}=\left\{\begin{array}{l}0 \\[3pt] 0\end{array}\right\} \end{eqnarray}
Setting the determinant of Equation (16) to zero, we obtain quadratic Equation (17) in terms of p, which is a reduced eigenvalue defined in dimensionless form by Equations (18) and (19).
\begin{equation}A p^{4}+B p^{2}+C=0, \end{equation}
\begin{equation} p_{1}=\frac{b v_{1}}{V}=\frac{b}{V}\left(\Gamma_{1} \pm i \Omega_{1}\right), \end{equation}
\begin{equation} p_{2}=\frac{b v_{2}}{V}=\frac{b}{V}\left(\Gamma_{2} \pm i \Omega_{2}\right), \end{equation}
where b is the semi-span and V is the velocity.
$\Gamma$
is the real part of the eigenvalue, which represents the modal damping, while the imaginary part
$\Omega$
represents the modal frequency. The advantage of using eigenvalues is that it allows an evaluation of the stability characteristics of the system in terms of frequency and damping. Table 1 presents the motion and stability characteristics of a system based on the real and imaginary parts of its eigenvalues.
Table 1 Stability characteristics based on eigenvalues(Reference Dowell14)
It is apparent that Equation (16) provides information only about the stability characteristics of the system at a given condition; it determines whether the system will be stable or unstable. However, it does not describe the stability boundary of the system. To yield the stability boundaries, an extended analysis based on this eigenvalue framework is carried out. The term
$B^2-4AC$
in Equation (17) determines the stability of the system, as it forms the imaginary part of p. Setting this term to zero results in Equation (20), which can be arranged in terms of the dynamic pressure q. The calculated dynamic pressure corresponds to flutter and can be used to obtain the flutter velocity
$V_f$
.
\begin{eqnarray}\left[m\left(K_{\alpha}-q S e C_{L \alpha}\right)+K_{h} I_{\alpha}-S_{\alpha} q S C_{L \alpha}\right]^{2}-4\left(m I_{\alpha}-S_{\alpha}^{2}\right)\left[K_{h}\left(K_{\alpha}-q S e C_{L \alpha}\right)\right]=0 \qquad\end{eqnarray}
3.0 FINITE ELEMENT MODELLING
The wing structure under study includes ribs made of Al 2024 alloy. Two ‘C-type’ spars run along the span. The outer half of each spar is made of a composite laminate, whereas the inner half towards the root is made of Al 7050 alloy. A thin sheet of composite laminate is wrapped around the wing to form the skin of the wing structure. An accurate finite element model (FEM) of the right wing was constructed (Fig. 4). The surface dimensions of all the components constituting the wing structure are significantly large compared with their thickness. To reduce the computational time and resources without compromising n the fidelity of the solution, the solid model can thus be accurately represented by a shell model. The thickness of the member is added as the thickness of the 2D shell elements. Isotropic properties are applied in the model for the metal structures, whereas orthotropic properties are assigned to the composite parts.
Figure 4. The reduced shell model generated for the wing.
The solution set-up and meshing were carried out in FEMAP, a modeller in NX Nastran developed by Siemens. A high-quality, quad-dominant structured mesh was created by manually defining equal numbers of elements on neighbouring edges (Fig. 5). A grid independence test with six different mesh sizes was carried out to select the optimum mesh size based on normalised frequency. As shown in Fig. 6, mesh convergence was obtained at 3,430,000 elements, which was thus chosen as the mesh size for all further analysis.
Figure 5. Finely generated mesh on ribs and skin cover of wing structure.
Figure 6. Grid independence is achieved at mesh size of approximately 3.4m elements.
4.0 MODAL ANALYSIS
Modal analysis is carried out to evaluate the dynamic response of the system in the frequency domain. It yields the mode shapes as a function of material properties, governing equations and boundary conditions. The frequency associated with this response is known as the modal or natural frequency. For simple structures, this can be done manually, but for complex structures, finite element analysis (FEA) is used, which is itself an intensive activity. Various methods are available to obtain the eigenvalue solutions in modal calculations, including the widely used determinant search, QR, QL, sectioning, subspace iterative and power (as well as inverse power) methods(Reference Bostic33). These methods can solve the complete set of equations describing the eigenvalue problem, in contrast to other, less rigorous procedures such as the modified power Lanczos method, which uses a recursive algorithm to find the desired set of eigenvalues. This method works by transformation of the larger and more complex eigenvalue problem into smaller and simpler trigonal problems, from which the eigenvalues and vectors of interest can be easily extracted. The general form of the eigenvalue expression for free vibration is given as Equation (21), where K, M,
$\omega$
and x follow conventional usage.
\begin{eqnarray}[K]\left\{x_{i}\right\}=\omega_{i}^{2}[{\rm M}]\left\{x_{i}\right\}\end{eqnarray}
The Lanczos solution naturally yields the highest values, while we are interested in the lowest few eigenvalues. To determine these, a shifting function
$\sigma$
is defined. Equation (21) is thus transformed into a reduced eigenvalue problem using a recursive technique.
\begin{equation}\omega^{2}=\frac{1}{\lambda}+\sigma,\end{equation}
\begin{equation}\lambda_{i}[K-\sigma {\rm M}]\left\{x_{i}\right\}=[{\rm M}]\left\{x_{i}\right\}, \end{equation}
\begin{equation}\left[A\right]\left\{x_{i}\right\}=\lambda_{i}\left\{x_{i}\right\},\end{equation}
\begin{equation}\left[A\right]=[K-\sigma {\rm M}]^{-1}[{\rm M}],\end{equation}
\begin{equation}\left[T_{j}\right]\left\{q_{i}\right\}=\lambda_{i}\left\{q_{i}\right\}\end{equation}
This method retains the previously computed results of the eigenvalue problem, which contributes to its high efficiency and rapid convergence rate. The result of the modal analysis for the first ten modes is shown in Figs. 7 and 8. The modal analysis feature integrated into the FE package was used to calculate these modes.
Figure 7. First bending mode (left) and first torsion mode (right) of FE model obtained by modal analysis in NX Nastran.
Table 2 Real and imaginary parts of eigenvalues of the studied UAV at different altitudes
Figure 8. Real and imaginary parts of first ten modes of vibration of FE model.
The negative real part of the complex eigenvalues represents the structural damping of the system. The contribution of the lowest modes of vibration is greatest, gradually decreasing for higher modes. Although the first ten modes of vibration are shown, the first six can actually be used for further flutter analysis. From the table, the most significant modes of vibration are the first and second bending, first torsion and first swaying mode.
5.0 RESULTS AND DISCUSSION
5.1 Eigenvalue analysis
All the parameters required for the eigenvalue framework described by Equations (16) and (20) are acquired from the FE model of the wing. The results obtained from the eigenvalue analysis using the analytical p-method are presented in Table 2. A positive real part of an eigenvalue indicates divergent behaviour, whereas a zero or negative value indicates a stable response. A non-zero imaginary part indicates an oscillatory response. From the eigenvalues analysis based on the p-method, we obtain a zero real part and non-zero imaginary part for all roots. This indicates that the wing experiences a stable oscillatory response at altitudes from 0 to 30,000 ft. From these results, it can be suggested that flutter does not occur in the designed operating conditions.
5.2 Extended eigenvalue analysis
The results of the eigenvalue analysis above suggest that the wing is free from flutter up to the speed limits defined by MIL STD 8870. However, the eigenvalue analysis does not provide the exact speeds at which flutter would actually occur. Therefore, an extended eigenvalue problem is solved to determine the flutter speeds. For this purpose, the parameters obtained from the computational model are inserted into Equation (20) to determine the flutter boundaries at various altitudes as shown in Fig. 9. The flutter speeds are found to increase with increasing altitude because of the lower density at higher altitudes, which reduces the magnitude of the aerodynamic forces acting on the wing.
Figure 9. The flutter speed versus altitude obtained from extended eigenvalue analysis. The design speed of the wing and the MIL STD speed are found to be well below its flutter speed at all altitudes.
5.3 Numerical validation (k-method)
The k-method introduces structural damping into the equations of motion, where the dissipative nature of the damping terms is defined as shown in Equations (29) and (30).
\begin{equation}m \ddot{h}+K_{h} h+S_{\alpha} \ddot{\alpha}=-L+D_{h}, \end{equation}
\begin{equation}S_{\alpha} \ddot{h}+I_{\alpha} \ddot{\alpha}+K_{\alpha} \alpha={\rm M}_{y}+D_{\alpha}, \end{equation}
\begin{equation}D_{h}=\bar{D}_{h} e^{i \omega t}=-i g_{h} m \omega_{h}^{2} \bar{h} e^{i \omega t}, \end{equation}
\begin{equation}D_{\alpha}=\bar{D}_{\alpha} e^{i \omega t}=-i g_{\alpha} I_{\alpha} \omega_{\alpha}^{2} \bar{\alpha} e^{i \omega t} \end{equation}
The structural damping factors
$g_h$
and
$g_\alpha$
depend on the structural and material properties, and these terms are assigned values in the range of 0.01–0.05 in the absence of any experimental data(Reference Hodges and Pierce13). Assuming a harmonic function, the final form of the equation to be solved using the k-method becomes
\begin{eqnarray}\left[-{\rm M} \omega^{2}+i B \omega+(1+i g) k-q Q_{{\rm M}, k}\right]\{U\}=0, \end{eqnarray}
where U is the displacement matrix comprising h and
$\alpha$
, the matrix B contains the aerodynamic state relations, and
$Q_{m,k}$
is a function of the reduced frequency (k) and Mach number (M). The k-method calculates the roots of Equation (16) with the Mach number, reduced frequencies and density as user-defined inputs. The k-method was run in NX Nastran software for the FE model as for the previous analysis, and flutter analysis was carried out at different altitudes ranging from sea level to 30,000 ft (in steps of 5,000 ft). The results of this analysis are shown in V–g and V–f diagrams in Fig. 10.
Figure 10. Frequency and damping of first six modes of flutter obtained by k-method, shown in V–g (top) and V–f (bottom) diagrams.
From the V–g diagram, it can be seen that the damping of the third, fourth, fifth and sixth modes changes sign from negative to positive, i.e. shifts from the stable, negative region to the unstable, positive region. For the given range of reduced frequency, the third mode shifts to the unstable region earliest. In the V–f diagram, corresponding to the point where the V–g trend changes direction, the frequency of the third mode starts to decrease continuously with increasing velocity. For better visualisation, the graphs for the first and third modes are shown in Fig. 11. In the V–g diagram, at
$V_a$
it can be seen that mode 3 starts to change direction and at
$V_f$
(20% higher speed than
$V_a$
), it crosses the axis dividing the positive and negative region. Correspondingly, in the V–f diagram, the frequency of the third mode remains constant up to
$V_a$
, after which it starts to decrease steeply and coalesces with the frequency of the first mode, which tends to remain unchanged with increasing speed. The diagram can be extended up to the point where the first bending and first torsional mode merge to give a flutter response, in which the dominant mode is torsional. Figure 11 shows the results of the analysis at sea level, whereas the results obtained at different altitudes are shown in Fig. 12.
Figure 11. V–f (left) and V–g (right) diagrams for the first torsion and bending modes. Teh frequency of the first and third modes approach each other with increasing velocity. In the V–g diagram, the damping of the third mode can be seen to cross the stability boundary.
Figure 12. Flutter speeds calculated at various altitudes using the p-, k- and pk-methods, revealing an increase with altitude.
5.4 Numerical validation (pk-method)
To validate the results obtained using the analytical p-method and computational k-method, the pk-method was also applied to carry out flutter analysis of the wing. While the p-method has established itself as the most accurate, the k-method is regarded as the fastest method(Reference Hodges and Pierce13). The pk-method is a combination of the k- and p-method, combining the benefits of both. It introduces structural damping into the p-method using a rate of decay
$\gamma$
.
\begin{equation} p=\gamma k+i k, \end{equation}
\begin{equation}\gamma=\frac{1}{2 \pi} \ln \left(\frac{a_{m+1}}{a_{m}}\right), \end{equation}
where
$a_{m+1}$
and
$a_{m}$
are the amplitudes of successive cycles. The final form of the equation becomes
\begin{eqnarray}\left[{\rm M} p^{2}+\left(B-\frac{\rho \bar{c}}{4 k} Q_{{\rm M}, V}\right) p+\left(k-q Q_{{\rm M}, V}\right)\right]\{u\}=0 \end{eqnarray}
Equation (34) is similar to Equation (31) of the k-method, except Q here is a function of M and velocity rather than of reduced frequency. The roots of Equation (34) are calculated for user-defined inputs of the Mach number and density for a range of altitudes. This analysis is run using Nx Nastran at different altitudes ranging from 0 to 30,000 ft in steps of 5,000 ft. A comparison of the critical flutter speeds calculated using all three methods, viz. the p-, k- and pk-methods, is shown in Fig. 12, revealing reasonably close agreement between all the methods. The k-method offers the most conservative results with the lowest flutter speeds at all altitudes. The p-method gives lower flutter speeds as compared with the pk-method up to 15,000 ft, after which the trend reverses and the pk-method gives lower speeds up to 30,000 ft. All the methods follow the same trend, in which the flutter speeds are found to increase with increasing altitude.
6.0 PARAMETRIC ANALYSIS
The results of these aeroelastic stability analyses of the wing using the three different methods suggest that the wing does not exhibit flutter at its design operating conditions. It is sometimes essential to further determine any parameters that may enhance the aeroelastic performance of the wing without significant design changes or weight penalties. While many methods of suppression are available, the purpose here is not to suppress flutter but to determine whether incorporating any design changes would improve the overall aeroelastic characteristics of the wing. Such parametric analysis can be done by altering the geometric parameters. The thickness of the spars and ribs is varied, and the corresponding flutter speeds are calculated. Increasing the thickness of the structural members implies more weight (Fig. 13). The changes in the gradient of the normalised flutter speeds at a particular altitude with respect to the thickness factors are shown in Fig. 14.
Figure 13. Normalized mass variation with thickness factor.
Figure 14. Normalised critical flutter speed versus thickness of structural members for rib, spars and ribs + spars.
The parametric study was carried out in three steps. First, only the thickness of the spars was changed while keeping all the other thicknesses constant, and the effects on the flutter speeds were calculated. In a second step, the thickness of ribs was changed while keeping the spar thickness constant. In the third step, the thickness of both the ribs and spars was changed, then the results from all three cases were compared.
Figure 15. Gradient of normalised flutter speed versus thickness factor for rib, spars and ribs + spars.
As seen from Fig. 13, changing the thickness of the spars results in the most significant change in the normalised mass factor, which means that thicker spars will take a higher toll on the weight of the overall structure. This is because the dimensions of the spars are greater than those of the ribs. It is evident that increasing the spar thickness increases the flutter speed linearly. However, changing the thickness of the ribs does not offer any advantages in terms of flutter speeds. Figure 15 shows a plot of the mass-normalised flutter speeds with respect to the thickness of the members for ribs only, spars only and the rib–spar combination. Clearly, increasing the thickness of only the ribs has no effect on the flutter speed gradient, which only contributes negatively for higher thickness values. The speed gradient for the rib–spar combination shows a better response, following the same trend as the case for spars alone. The spar-alone configuration thus offers the most benefit, with a positive gradient even for members with lower thickness than in the actual design. Wings spars are designed to take bending as well as torsional loads, because the stiffness of the spar offer better resistance against deformation than the ribs. Therefore, based on the results of this parametric study, it is only reasonable to state that the wing flutter speeds are significantly influenced by the spar thickness. The aircraft wing can be modelled as a cantilever beam with a stiffness given as a function of the Young’s modulus (E), Shear modulus (G), polar moment of inertia (J), second moment of area (I) and the boundary conditions.
\begin{eqnarray}{\rm K} \propto {\rm E}, {\rm G, I, J} \end{eqnarray}
J and I depend on the geometric distribution of the beam cross section. These inertial and moment properties of the wing vary when changing the structural thickness of the spars, thus having a direct impact on the overall stiffness and hence the flutter speeds of the wing.
7.0 CONCLUSION
This study set out to investigate a novel approach to determine the aeroelastic characteristics of a modern HARW design for UAVs. The study used parameters interchangeably between the analytical and computational analyses. The k-method proved to be most conservative, as it gave the lowest flutter speeds, while the pk-method gave the highest flutter speeds at altitudes from sea level to 15,000 ft, and the critical flutter speeds obtained using the p-method were highest from 15,000 to 30,000 ft. All three methods appear to follow the same trend. The Lanczos method coupled with a recursive technique was applied to obtain the first ten significant modes of vibration. The extended eigenvalue solution was necessary to obtain the flutter boundaries, with the flutter speed as the objective function. It was found that the wing will exhibit flutter due to the coalesce of the first bending and torsion mode. As per the requirements of MIL STD 8870, the wing should be free of flutter at speeds beyond 15% of the maximum design flight speed. This research was extended to include a sensitivity analysis in which the thicknesses of the spars and ribs were changed and the corresponding flutter speeds determined. Changing the thickness of the ribs offered no benefit, as no improvements in the flutter response were observed; instead, it contributed adversely at higher thickness factors. On the other hand, increasing the thickness of the spars translated into higher flutter speeds. These differences can be explained in part by the better resistance of thicker spars to bending and torsional loads. The results of the parametric study indicate that the aeroelastic response of the wing is more sensitive to the bending and torsional stiffness of the spars.
While this analysis is carried out with reasonable fidelity and its validation is also presented herein, some factors that were assumed or neglected in this study could be considered in future analyses. A separate study to investigate the differences obtained in the results when incorporating such factors could help to justify their exclusion or, if required, suggest methods to include them if their effects are significant enough to be considered in future analyses. For example, damping parameters from GVT could be used instead of historical values. The basic wing analysis is carried out here without considering the effects of control surfaces. The flutter response of a wing at different speeds and altitudes with the movement of control surfaces could thus also be studied. The stiffness of the wing–fuselage attachment is also neglected, as a wing-alone configuration is studied with the wing modelled as a cantilevered beam. Moreover, analysis of the entire UAV geometry would capture the flutter and buffeting responses of the horizontal tail as well, which are significantly affected by the downwash from the wing. Nonlinear effects are not taken into account in this analysis but may arise from stall flutter and the free play of control systems. Non-linearities in the structure and aerodynamics could also be taken into account as part of future studies.
