1. Introduction
Aeroelasticity is a multidisciplinary field consisting of aerodynamics, structural mechanics, dynamics, and their interactions. By assuming linear aerodynamics and structural dynamics, classical theories in aeroelasticity are able to predict dynamic aeroelastic instability, i.e., coupled flutter, which may define the stable region of an aircraft flight envelope. These linear theories can be implemented in the frequency or time domains. As they are extremely efficient numerically, frequency domain solutions are generally preferred. However, they contain one important drawback, in particular, in the context of unsteady aerodynamics (Edwards and Wieseman [
1]). To consider unsteady aerodynamics effects, the generalized aerodynamic forces (GAF) are formulated in the frequency domain. The GAF, symbolized by
Q(
ik), are tabulated for a number of sets of discrete reduced frequencies,
where
i is the unit imaginary number,
is the characteristic chord length and
V is the airspeed. The GAF is not a continuous analytical expression that could be imbedded into the structural force terms to form the full aeroelastic equations, leading to the requirement of an iterative procedure such as frequency matching (Wright and Cooper [
2]) because the value of the reduced frequency of
k as input to the flutter equation is part of the flutter solution.
Doublet-lattice method (DLM) is the most commonly used method for flutter analysis and certification in industrial practice. DLM is based on the compressible acceleration potential theory for thin wing geometry and it cannot account for the wing thickness nor capture recompression shocks or boundary layer separations. However, flutter and other aeroelastic problems are often observed in jet aircraft flights, where nonlinear aerodynamic characteristics are present due to local supersonic flow regions terminated by shocks. Since computational fluid dynamics (CFD) becomes mature in capturing shock waves and simulating shock wave–boundary layer interactions, CFD is increasingly used for computational aeroelasticity. In order to achieve this, a number of researchers coupled CFD with computational structural dynamics (CSD) [
3,
4,
5]. The coupled CFD-CSD approaches effectively showed capability in accounting for the nonlinear aeroelastic effects, including the transonic dip due to compressibility effects and limit cycle oscillations (LCOs). In the fully coupled CFD-CSD methods, the simulations are carried out in the time domain and the structural displacement responds instantly to the forces exerted by the fluid, which requires large computational cost.
As an efficient industrial practice, CFD is also used to correct or replace the DLM in flutter analysis, forming a hybrid linear–nonlinear or time–frequency domain methods. Kaiser et al. applied a CFD-corrected DLM to account for nonlinear effects in a gust-encounter investigation [
6]. Fleischer and Breitsamter used CFD methods for the computation of GAFs [
7]. A linear frequency domain solver was used for harmonic small-disturbance solutions in [
6] while a time-marching Euler method was used to generate the GAFs in [
7]. In the latter case, the CFD results obtained in the time domain were transformed to the frequency domain using Fourier analysis for calculation of the GAFs used for the following flutter analysis.
The aforementioned fully coupled CFD-CSD approaches are used in the time domain while the hybrid linear–nonlinear methods use CFD to correct or replace the DLM and its corresponding aeroelastic stability analysis for complex eigensolutions still remains in the frequency domain, the same as in the linear flutter analysis methods. In the widely used linear methods, there are three main steps: 1. Structural modal analysis to generate the real eigenvectors as the vibrational mode shapes; 2. Aerodynamic modelling to generate the aerodynamic influence coefficients (AICs) and GAFs; 3. Aeroelastic stability analysis searching for complex eigensolutions to determine flutter speed. The aforementioned hybrid methods using CFD improve the aerodynamic modelling in Step 2. The present study is focused on numerical stabilizations of the aeroelastic stability analysis process in Step 3.
There exist a number of flutter analysis approaches to solve the frequency-matching problem. Among them, the K-method and the PK-method are most commonly applied. Both methods are based on simple harmonic motion assumptions, which are acceptable at the flutter boundary. However, they are not valid below and above the flutter boundary. Consequently, the methods may predict the same flutter speed and frequency, but with an inaccurate subcritical behavior. It is known that damping predicted using the K-method is not reliable, c.f., Section 5.4.2 in Hodges and Pierce [
8]. Therefore, the PK-method is often selected because of its reliability.
The PK-method has wide applications, which are exemplified by the following two recent test cases. Dinulović et al. [
9] applied the PK-method to flutter analysis of tapered composite fins based on panel vortex methods. Rho et al. [
10] performed a hybrid frequency and time domain analysis to predict aeroelastic phenomena such as flutter, divergence and buffet of a hammerhead launch vehicle, where CFD was employed to calculate the unsteady aerodynamic forces for the vehicle in harmonic vibrating motions in the time domain. The aerodynamic forces were then transferred to the frequency domain for flutter analyses using
v-g and PK-methods.
In the PK-method, the Mach number, altitude and airspeed are fixed for each iteration, while the reduced frequency is set with an initial guess that is normally a low value. Starting from the mode with the lowest structural frequency, iterations are performed with the reduced frequency of the aerodynamics for each mode until the reduced frequency converges. Then, the iteration process is repeated for the next mode at a higher frequency until a converged solution is obtained for all the structural modes selected by the user.
Since an aeroelastic stability analysis of a real-world complex aircraft configuration may involve up to 100 vibrational modes, severe mode switching is often observed in the flutter analyses when using the commonly applied PK-method of the commercial code NASTRAN. Severe misleading mode switching may cause misidentification of aeroelastic phenomena. As a result, factors leading to instabilities such as flutter and divergence may be improperly interpreted and incorrectly understood, and attempts to make engineering decisions (such as stores clearance by flutter similarity) to improve flight performance could be misguided and likely unsuccessful.
To alleviate the severity of the aforementioned mode switching, mode-tracking methods based on eigenvalues and eigenvectors as a post-processing module have established a more reliable visualization of the mode behavior in flutter analysis. Nevertheless, issues with mode switching have not yet been fully resolved. We believe that mode switching is inherently correlated to the PK iteration process. As pointed out by Eldred et al. [
11], the PK-method has a feature where there may be more modes for which the aerodynamics match the complex eigenvalue (root) than the pre-defined structural ones at a given airspeed. The extraneous eigenvalues cause problems for the mode-tracking process since the algorithms may mistakenly lock onto one of these roots if it is similar to the true one.
This is mainly attributed to the fact that in addition to the aeroelastic eigenvalues that emerge from the structural natural frequencies, there can be aerodynamic eigenvalues that emerge from the aerodynamic model itself. Dashcund [
12] developed both theoretical and experimental models to study active suppression of flutter of a classical bending–torsion wing by adding a control surface. By comparing the predicted results to the experimental results, it was determined that some aeroelastic roots resulting in divergence were from neither structure nor control system. They were derived from the aerodynamic model instead. By examining the eigenvectors of a two-dimensional (2D) wing that had rigid airfoil sections with permit for rotation only, Heeg [
13] reported that the least stable real mode originated in the aerodynamics and became the source of static divergence. This finding was proven by their wind-tunnel experiments using a quasi-2D wing. More recently, Vedeneev [
14] confirmed that when applying unsteady aerodynamics to a two-degree-of-freedom (two-DOF) bending–torsion wing, the structural eigenfrequencies became damped while the divergence mode separated from a continuous spectrum that existed in the aeroelastic system due to the wake behind the wing that was embodied in the unsteady aerodynamic model.
There exist a number of mode-tracking techniques. Desmarais and Bennett [
15] relied on the shape of the characteristic polynomial for their
v-g analysis and used Laguerre iteration to converge from a previous to a corresponding current eigenvalue. van Zyl [
16] correlated the modes based on complex inner products between current and previous right eigenvectors. Eldred et al. [
11] proposed complex higher-order eigenpair perturbations (C-HOEP) and complex cross-orthogonality check (C-CORC) methods. However, they concluded that in their PK-method analysis, van Zyl’s method outperformed the C-CORC and C-HOEP methods. Recently, Ren [
17] developed a predictor-corrector scheme based on the perturbation theory of eigenvalues for mode tracking and applied it to the piecewise quadratic interpolation (PQI) method proposed by Goodman [
18]. The correction scheme adaptively adjusted the speed interval between the current and next speed steps by comparing the estimated eigenvalues and sorted eigenvalues. Similarly, Quero et al. [
19] developed a mode tracking algorithm by adapting the interval of the parameter of interest (i.e., airspeed in PK-method) and applied it to the p-L flutter analysis method generalized by Quero et al. [
20]. Nevertheless, simply using smaller airspeed intervals in the PK-method did not improve the results for the real-world complex configurations discussed in the present study.
Moreover, solutions with eigenvalues (roots) less than the vibratory mode numbers were also observed in the present study. When the numbers of modes and roots do not match, the post-processing algorithms for mode tracking are no longer effective and can be meaningless. Thus, investigations on the PK iteration process are needed to improve the numerical stability of the flutter analysis methods instead of pure re-ordering (mode tracking). As expected, this study confirmed that the initial guess and the iteration procedure of the eigensolution are crucial to obtain stable and accurate flutter trends.
It should be noted that a promising
g-method was proposed by Chen for reliable damping prediction [
21]. Gu and Yang further reformulated and solved the damping iterations in a modified PK-method [
22]. The methods were well illustrated for wings with less than 10 vibrational modes. However, it has not been confirmed whether the
g-methods can improve the issues with severe mode switching for real-world complex applications involving dozens of structural modes.
To make use of the structural and aerodynamic data outputted from the commercial finite element software NASTRAN as input to our in-house flutter solver FLUTQ, a communication interface between NASTRAN and FLUTQ was developed in this study. Stabilization techniques were developed and implemented for the PK-method, including mode tracking between airspeeds, mode matching to lock eigenroots onto the aerodynamics, a hybrid scheme for the initial guess of the reduced frequency k, a deferred scheme used for the PK-iteration procedure, and a modified g-method using damping iterations. Consequently, the applicability of the in-house solver FLUTQ was extended. As a result, the numerically improved aeroelastic stability analysis effectively minimized the severity of the misleading mode switching often observed in flutter analyses for real-world aircraft. The more reliable prediction and higher clarity of the flutter modes led to less time and effort to understand and interpret the flutter analysis results, with increased confidence that the predictions were valid for reducing risks in flight.
It should be mentioned that one can extract true aeroelastic eigenvalues by allowing the reduced frequency in the unsteady aerodynamic models to be complex. Edwards et al. [
23] employed the generalized Theodorsen function for the aerodynamic loading of a 2D thin wing with three-DOF motions. Instead of the assumption of simple harmonic motion, the transient response was solved and the true aeroelastic eigenvalues including the damping ratios and frequencies were obtained. Liska and Dowell [
24] applied generalized or complex Theodorsen function in an unsteady aerodynamic model to the aeroelastic equations of a two-segment folding wing. The complex Theodorsen function depends on the motion’s frequency and damping, and describes the magnitude and phase of circulatory effects in unsteady aerodynamics. Although exact eigensolutions were obtained in these studies, these methods are computationally expensive. Moreover, the use of analytical aerodynamic models is neither feasible nor realistic for applications to real-world complex configurations.
It should be noted that an initial version of this work was presented in [
25]. The major difference between the initial version and this manuscript is that both the eigenvector-based mode tracking technique and a modified g-method were recently implemented in the in-house solver FLUTQ. The FLUTQ results presented in the initial version were neither mode-tracked using the eigenvector-based approach nor obtained with a damping iteration using the modified g-method. In next sections, the numerical methodologies will be briefly described first, followed by demonstrations of the improvement to the results.
4. Demonstration and Discussion of the Computed Results
The developed methodologies were first verified against the NASTRAN tutorial example HA145B introduced by Rodden and Johnson [
27]. This example presents a flutter analysis of a BAH jet transport wing. The structural, geometric, and aerodynamic parameters are provided by the tutorial in the form of NASTRAN input.
Figure 4 shows the BAH wing planform extracted from the NASTRAN user’s guide. The aerodynamic model was divided into six partitions along the semi-span and four boxes chordwise. No motion was assumed at the fuselage centre. Eleven structural grid points are located on the wing, resulting in 10 coupled frequencies of the cantilever wing. The wing properties are reported in detail by Rodden and Johnson [
27]. The wing was assumed to be flying in incompressible air (Mach 0.0) and at a dynamic pressure
q = 4.0075 psi. The subsonic DLM in NASTRAN was used for the flutter analyses by the PK- or
g-methods. The NASTRAN tutorial input files were modified in this study to output the aerodynamic matrix to allow application of the in-house FLUTQ solver.
The developed methodologies were applied with the 10 vibrational modes specified. The modal damping and frequency plots obtained with three different techniques are shown in
Figure 5, where the top diagrams were from the NASTRAN PK-method, the middle diagrams were from the FLUTQ PK-method, and the bottom diagrams were from the FLUTQ
g-method. All three sets of results suggested that a lowest flutter speed is a coupled first-bending first-torsion mode flutter and found in
Vf ≈ 1054 ft/s and
Ff = 3.09 Hz at
g(G) = 0. The
v-f and
v-g trends were clear with no appearance of mode switching, indicating that all the methods delivered nearly identical flutter results for this simple configuration with only a limited amount of vibrational modes considered. However, the following test cases will exhibit a much more challenging situation for real-world complex aircraft configurations.
After verification using the simple wing model, the improved methodologies were applied to four real-world jet aircraft configurations with underwing and possible tip-mounted stores: two symmetric and two asymmetric configurations as illustrated in
Figure 6. The stores represent any objects such as external fuel tanks, weapons, and special electronic sensors or devices. The results were compared with those obtained from NASTRAN.
The flutter analysis model of the real-world jet aircraft was a linear model. A whole-aircraft model was used, which allowed flutter analysis for asymmetric configurations and better prediction of anti-symmetric flutter phenomena. The structural model was a flexible beam or stick model where the fuselage, wing, wing controls and wing tip-mounted stores were modelled by flexible beams and concentrated mass elements. Wing pylons and stores were modelled by a small number of discrete scalar stiffness and concentrated mass elements. The rigid empennage was simply modelled as additional fuselage point masses with no flexibility.
The aerodynamics of the aircraft were modelled using the DLM. The one-side single wing consisted of 140 panels. The fuselage was modelled by five elliptical-section slender body elements and five cylindrical interference elements. The empennage and underwing stores were not modelled aerodynamically. The tip-mounted store was modelled as an open-ended octagonal cylinder with 184 panels. The flutter analyses assumed the following:
- (1)
Mach 0.95 aerodynamics;
- (2)
matched altitude flutter solutions, i.e., the air density varied at each speed to match the assumed Mach number;
- (3)
80 flutter solution modes for the whole aircraft;
- (4)
no structural damping;
- (5)
calculation of the g = 0.0 flutter crossing.
The first assumption with Mach 0.95 aerodynamics is obviously rough. As NASTRAN relies on the subsonic DLM, it is not fully possible to determine the unsteady aerodynamic loads to a sufficient accuracy level for transonic conditions. In a study on control surface aerodynamics, Roughen et al. [
37] confirmed that DLM results for phase angle agreed well with both CFD and wind tunnel data except for in supersonic regions and cases with high reduced frequency. DLM predicted too small a phase lag in supersonic regions as expected because linear potential theory is not valid in regions where compressibility has a significant influence. Despite the possible inaccuracy, DLM as a common industrial practice has a long history. Triplett used DLM in aerodynamic modelling for F/A-18 wing-store tip-missile flutter at Mach 0.9 and 0.95 [
38]. Although the current study is focused on the numerical procedure stabilization in searching for complex eigensolutions, as discussed in the introduction, CFD corrections may improve the accuracy, which should be considered in future studies.
4.1. Effectiveness of the Mode-Tracking Algorithms
Figure 7 shows impacts of the van Zyl [
16] method when it was applied as a post-processing mode-tracking procedure to the Symmetric 1 configuration of a real-world complex jet aircraft. The PK-method flutter analysis was conducted using 80 vibrational modes. The first six modes were rigid-body modes. For clarity, the figure displays only the first 14 flexible modes (modes 7–20). The airspeed was arbitrarily normalized in the velocity–damping (
v-g) and velocity–frequency (
v-f) plots. The computed damping coefficient G(
g) indicates the artificial structural damping required to damp the oscillation caused by the interaction between the aerodynamics, structural mechanics, and dynamics. Negative values indicate a stable aeroelastic state while positive values indicate flutter occurrence that requires structural damping to suppress the oscillations. The predicted lowest flutter speeds were
Vf = 77.35% from NASTRAN and
Vf = 77.54% using FLUTQ, with a frequency of 5.69 Hz at
g = 0. The principal flutter modes were mode 12 (Antisymmetric Wing First Bending) and mode 14 (Antisymmetric Underwing-Store Pitch). As shown in the figure, the mode-tracking algorithm as a post-processing procedure did not change the flutter prediction results, which was expected. However, it sorted the results and significantly reduced the amount of severe mode switching, thus alleviating confusion in the interpretation of the flutter results. Note that the FLUTQ PK-method results showed less mode switching when compared with the NASTRAN PK-method results.
The mode-tracking mechanisms alleviated the severity of the misleading mode switching for some conditions. However, they do not work properly for all complex configurations.
Figure 8 demonstrates the ineffectiveness of the eigenvector-based van Zyl mode-tracking algorithm when it was applied to the NASTRAN PK-method results for an asymmetric configuration of the same parent jet aircraft. Although the mode-tracking algorithms reduced mode switching to some extent, the issue related to the mode switching was still obvious. A closer look at the mode-tracking process revealed that the ineffectiveness of the mode-tracking algorithms was mainly attributed to the inaccurate eigensolution. As marked by the red circle in
Figure 8b, the solution of a mode was missing in the original eigensolution at airspeed
V = 72.5%. The missing solution was inherent to the PK iterations, and the post-processing mode tracking was unable to correct it. As confirmed by
Figure 8c, this problem was solved effectively by using the in-house flutter analysis solver. For the Asymmetric 2 configuration, the predicted lowest flutter speeds were
Vf = 61.30% from NASTRAN and
Vf = 61.33% using FLUTQ, with a frequency of 8.08 Hz at
g = 0. The principal flutter modes were mode 14 (Antisymmetric Fuselage First Bending) and mode 15 (Right Hand Underwing-Store Pitch).
4.2. Solution Stabilization for the PK-Iteration
Figure 7 and
Figure 8 also demonstrate the effectiveness of the solution stabilization applied to the PK-method using the in-house solver FLUTQ. The structural and aerodynamics matrices used for the FLUTQ flutter analyses were the same as those used for the respective NASTRAN solutions. An interface was programed to extract these data from NASTRAN to be employed for the flutter analyses using the FLUTQ PK-method. As shown in the figures, the FLUTQ solutions were extensively improved, and no severe misleading mode switching appeared until the post-flutter region.
Figure 9 further demonstrates the effects of the individual numerical techniques. For simplicity, the figure only displays the first 10 flexible modes. The results shown in
Figure 9a demonstrate a test case for which no mode matching-and-locking procedure was used in the PK iterations, confirming that mode matching-and-locking is a key element to stabilize the solution.
The effects of the deferred scheme and the hybrid initial solution on the results were not as obvious as the mode matching-and-locking mechanism. In general, Equations (16) and (13) worked well for most test cases (when other techniques developed in this study were applied). However, introducing the deferred correction scheme in Equation (17) and the hybrid initial solution approach in Equation (14) using under-relaxation parameters alleviated mode switching in the post-flutter region, as evidenced in
Figure 10. Mode 15 appeared as a hump flutter mode that crossed
g = 0 shallowly and then re-crossed to the stable side of the
v-g plot. Hump flutter modes are most commonly observed in store flutter. The damping of mode 15 was far below 0.02 (an assumed inherent structural damping commonly used in industry) and was believed unlikely to produce explosive flutter. Therefore, the predicted flutter speed was
Vf = 75.54% with a frequency of 7.29 Hz when crossing
g = 0. The principal flutter modes were mode 12 (Antisymmetric Underwing-Store Pitch) and mode 14 (Antisymmetric Fuselage First Bending).
Figure 10a,b show more mode switching, which was mainly attributed to numerical divergences in the post-flutter region when compared to
Figure 10c. The under-relaxation parameters provided an additional option to control the numerical stability in the PK iterations.
Note that except for the one under investigation, all other numerical stabilization techniques developed in this study were applied in the demonstration tests of the effects of the individual numerical techniques, which is slightly different from the work presented in our initial publication [
25]. Consequently, the results presented in this paper may have a better looking when compared with those presented in the previous paper.
4.3. g-Method Results
The implemented
g-method was first tested with a NASTRAN tutorial example using the HA145B wing introduced by Rodden and Johnson [
27]. After successful verification, the method was applied to flutter analyses of real-world jet aircraft configurations, as performed for the PK-method, to demonstrate its capability.
Figure 11 and
Figure 12 verify the
g-method implemented in the FLUTQ solver against the PK-method. The predicted flutter speeds for the Symmetric 2 configuration shown in
Figure 11 were
Vf = 75.53% using the
g-method vs.
Vf = 75.54% using the PK-method, with a frequency of 7.29 Hz when crossing
g = 0. The principal flutter modes were mode 12 (Antisymmetric Underwing-Store Pitch) and mode 14 (Antisymmetric Fuselage First Bending). The predicted flutter speeds for the Asymmetric 1 shown in
Figure 12 were
Vf = 82.34% with
Ff = 5.70 Hz using the
g-method vs.
Vf = 82.33% with
Ff = 5.71 Hz using the PK-method. The principal flutter modes were mode 12 (Left Hand Underwing-Store Pitch) and mode 14 (Antisymmetric Wing First Bending). For the results, both eigenvalue- and eigenvector-based mode-tracking mechanisms were applied to both the PK- and
g-method results. As shown in the figures, the implemented FLUTQ
g-method predicted reasonable flutter results. In general, the FLUTQ
g-method applied to the real-world jet aircraft slightly further reduced the severity of the misleading mode switching in the post-flutter region compared to that of the PK-method. This is because the aerodynamic damping is more realistic when the
g-method is used.
Parametric studies with different numerical settings were also performed for the numerical process in the g-method. In general, and as observed in the previous PK-method results, the combined eigenvalue- and eigenvector-based method performed better than either of the individual methods. When using the combined method, the eigenvalue-based mode tracking was performed first, followed by further sorting using the eigenvector-based method in the post-processing.
Slightly different from the PK-method, the deferred correction used for solution stabilization in the g-method did not show consistent performance, as observed in the previous PK-method results. In particular, the deferred correction technique improved the results only for some configurations. Therefore, it is recommended to switch off the deferred correction technique when the g-method is applied. This will save some computational time because the deferred correction using an under-relaxation parameter may require more iterations for convergence.
As discussed earlier, Equation (20) is valid for
g << 1; therefore, |2
γ| ≤ 0.02 was applied to
g =
γω when calculating the damping terms highlighted in Equation (23). To confirm the influence of this parameter, enlarged values of the parameter were studied. In general, |2
γ| ≤ 0.02 resulted in the clearest results over the ranges |2
γ| ≤ 0.04 and |2
γ| ≤ 0.1. Consequently, |2
γ| ≤ 0.02 or, in other words, |
g| ≤ 0.01
ω for Equation (32), is recommended to bound the damping terms in
g-method flutter analyses for real-world jet aircraft configurations. Physically, this value used for damping bounding is well aligned with the structural damping coefficient
g = 2
γ = 0.02 that is often used to calculate flutter crossing to determine the flutter speed for real-world aircraft configurations when considering realistic structural damping as discussed in
Section 3.3.2 when describing the methodologies.