Flutter analysis and mode tracking of aircraft model based on piecewise interpolation method

2.1. Piecewise quadratic interpolation of the generalized aerodynamic forces

The dynamic equation used for flutter analysis can be transformed into the following equation by using the Laplace transform [7]:

where $\stackrel{-}{\mathbf{M}}$, $\stackrel{-}{\mathbf{D}}$, and $\stackrel{-}{\mathbf{K}}$ are the generalized mass, damping and stiffness matrices, respectively. $L$, ${\rho }_{\mathrm{\infty }}$, and $V$ are the reference length of half-chord, air density and the flow speed, respectively. The generalized aerodynamic force ${\mathbf{Q}}_s\left(p\right)$ depends on the eigenvalue $p$, and can be computed via the doublet-lattice method [8]. The complex eigenvalue can be written as $p=g+ik$ where $g$ is the damping and $k$ is the reduced frequency.

Because ${\mathbf{Q}}_s\left(p\right)$ is a nonlinear asymmetric matrix, the direct solution of Eq. (1) can be computed by using the iterative methods. Typically, the eigenvalue problem is solved by introducing a parameter $\stackrel{-}{p}$ and computing the generalized aerodynamic force ${\mathbf{Q}}_s\left(\stackrel{-}{p}\right)$. Combing with a search or optimization procedure, the linear eigenvalues are computed and the eigenvalues equal to the parameter $\stackrel{-}{p}$ are saved as the solutions of the nonlinear Eq. (1). However, the iteration method not only increases the calculation time but also probably causes the problems with convergence. To improve the efficiency of the solution procedure, in this study, the relationship between ${\mathbf{Q}}_s\left(p\right)$ and $p$ is approximated by the following piecewise quadratic function:

where $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$ are the complex-valued coefficient matrices. $j$ and $b$ are the number of segments and the breakpoint of the reduced frequency range. The advantage of the piecewise quadratic function of ${\mathbf{Q}}_s\left(p\right)$ is that it converts the problem of nonlinear iteration procedure into a piecewise, linear quadratic eigenvalue problem. Moreover, because the Eq. (2) fulfills the Cauchy-Riemann equation [9] [defined in Eq. (3)], Eq. (2) can be constructed along the imaginary axis over the full complex plane. At the same time, it can be seen that Eq. (2) provides a second-order approximation of the variation with $g$:

The procedure of piecewise quadratic interpolation of aerodynamic forces can be illustrated as follows. At first, the generalized aerodynamic matrices $\mathit{Q}\left(k_i\right)$ are computed at different reduced frequencies $k_i$ ($i=\mathrm{}$1, 2,…, $n$) by using the doublet-lattice method. Then, the $n-$1 breakpoints $b_i$ for the piecewise interpolation are chosen as:

For each segment, three conditions must be satisfied. At first, the generalized aerodynamic forces on the boundaries must be matched exactly. Then, the slope values on the boundaries need to be matched. At last, the discrete values are interpolated exactly. The mathematical descriptions of the three conditions can be written as follows:

For the first and last segment, the conditions are modified to interpolate the generalized aerodynamic forces ${\mathit{Q}}_s\left(k\right)$ at the first and last reduced frequencies, respectively. For the first segment, the first condition [defined in Eq. (7)] is modified by the following formula:

For the last segment, the third condition [defined in Eq. (9)] replaced by:

Because the conditions are related between two neighboring segments, these conditions can be expressed in the matrix form by using the partitioning and assembly of matrices and the coefficient matrices ($A_j$, $B_j$, $C_j$) can be solved through the matrix transformation. For instance, the coefficient matrices ($A_j$, $B_j$, $C_j$) of Eq. (2) to be solved at the first and last segments can be written as:

where ($({\mathbf{A}}_{ts}{)}_j$, $({\mathbf{B}}_{ts}{)}_j$, $({\mathbf{C}}_{ts}{)}_j$) are the coefficient matrices for the $j$th segment, $t$ and $s$ are the row and column indexes in the matrix. ${\left(Q_{ts}\right)}_{ki}$ is the generalized aerodynamic force at $k_i$ reduced frequency.

2.2. Predictor-corrector based mode tracking

Once the piecewise quadratic function of ${\mathbf{Q}}_s\left(p\right)$ is available, Eq. (1) can be recast as:

Because ${\mathbf{Q}}_j\left(p\right)$ is only valid within the $j$th segment, the eigenvalues of Eq. (13) which yield the condition of $\mathrm{I}\mathrm{m}\left(p\right)\in \left[b_j,b_{j+1}\right]$ is kept. If there is no any eigenvalue satisfying the condition, it needs to continue to solve the eigenvalues by using the next segment function of ${\mathbf{Q}}_{j+1}\left(p\right)$ until the corresponding eigenvalues of each mode at a given flow speed are obtained. The process is implemented repeatedly at each flow speed and the flutter boundary and mechanism can be determined by drawing the $v$-$g$ and $v$-$f$ diagrams. It is noted that if the corresponding eigenvalues of all modes with $\mathrm{I}\mathrm{m}\left(p\right)\in \left[b_j,b_{j+1}\right]$ belong to different segments, the correspondence between the modes and eigenvalues due to the variation of flow speeds should be determined accurately. However, if there are two or even more eigenvalues kept in the same segment, then the mode tracking phenomenon will occur.

The issue of eigenvalue tracking can be handled by means of a Predictor-Corrector scheme which is based on the perturbation theory of eigenvalues. The forecast method is developed so as to track all modes during the flutter analysis. At the $m$-th flow speed $V_m$, the Predictor-Corrector scheme uses the eigenvalues ${\mathbf{p}}_i^m$ and the corresponding left and right eigenvectors to estimate eigenvalues $g^{m+1}$ at the next flow speed $V_{m+1}$. The estimated eigenvalues can be computed as follows:

where:

${\mathbf{\psi }}_i^m$ and ${\mathbf{\varphi }}_i^m$ are the left and right eigenvectors at the flow speed $V_m$, respectively. Once the estimated eigenvalues are available, the eigenvalues ${\mathbf{p}}^{m+1}$ can be sorted according to the estimated values. According to the minimal Euclidean metric between ${\mathbf{g}}^{m+1}$ and ${\mathbf{p}}^{m+1}$, the correspondence between the modes and the eigenvalues can be determined accurately. The minimal Euclidean metric can be computed as follows:

where:

In order to avoid the failure of the mode tracking scheme when applied the case with a larger step size of the flow speed, a threshold level should be set such as $\epsilon =\mathrm{}$0.001. For the large step size of the flow speed, when $d\left({\mathbf{g}}_i^{m+1},{\mathbf{p}}_i^{m+1}\right)>\epsilon$, it is needed to decrease the step via the factor $T$ ($T>1$), and recalculate the eigenvalues and sort them. The block diagram of the PQI method with predictor-corrector based mode tracking for flutter analysis is shown in Fig. 1.

3.1. Subsonic flutter analysis of a 15-degree sweptback wing

The 15-degree sweptback wing model is a standard example of flutter analysis in MSC Nastran. The modal parameters, such as the generalized mass, stiffness, and mode shapes were obtained from the MSC Nastran result file. The unsteady aerodynamic forces are computed by using the doublet-lattice method. Fig. 2 shows the configuration of the 15-degree sweptback wing model. For the flutter analysis of 15-degree sweptback wing model, the first four dominant modes were used. The mode shapes and natural frequencies of the four dominant modes are shown in Fig. 3.

Fig. 2Configuration of the 15-degree sweptback wing model

In the case study, the mode tracking problem and the accuracy of the PQI method with mode tracking are illustrated and verified The reduced frequency range is selected as 0.0-1.0 and divided into 50 segments to construct the piecewise quadratic function for each segment. The flow speed range for flutter analysis is from 10 m/s to 220 m/s and the step size of flow speed is 5 m/s. Fig. 4 shows the $v$-$g$ and $v$-$f$ diagrams computed via the PQI method without mode tracking. The computed dampings and frequencies of the second and third modes exchange each other all of a sudden. To clearly illustrate the reason for such a phenomenon, Table 1 gives the imaginary parts of the eigenvalues and the corresponding reduced frequency segments. As shown in Table 1, when the flow speed increases to 60 m/s, there are two roots with $\mathrm{I}\mathrm{m}\left(p\right)\in \left[\text{0.05,}\text{0.07}\right]$ kept in the same reduced frequency segment (0.05-0.07) but not sorted. This phenomenon can also been found when the flow speed increases to 80 m/s, 85 m/s, and 90 m/s, respectively.

Fig. 3Four dominant natural modes of the 15-degree sweptback wing model

a) The first mode (bending), 34.24 Hz

b) The second mode (torsion), 196.2 Hz

c) The third mode (torsion), 240.34 Hz

d) The fourth mode (bending), 563.01 Hz

Fig. 4The solutions of PQI method without mode tracking

a)$v$-$g$ plot

b)$v$-$f$ plot

Fig. 5 shows the $v$-$g$ and $v$-$f$ diagrams computed by the PQI method with mode tracking and the results are also compared with the P-K method. As shown in Fig. 4, the mode tracking problem of the PQI method is well handled by using mode tracking. Moreover, Fig. 5 also illustrates that the flutter boundary predicted via the PQI method gets a good agreement with the P-K method. However, there are small disagreements between the solutions of the damping ratios when the flow speed is lower than 185 m/s. When the flow speed is equal or greater than 185 m/s, there is a sudden jump of the damping of first mode computed via the P-K method. The reason is that the computed damping and frequency of the mode 1 are exchanged with the aerodynamic lag roots all of a sudden [10]. In contrast to the P-K method, the PQI method tracks all eigenvalues very well. The reason is that the PQI method provides a second-order approximation of the variation with $g$. It is noted that the eigenvalues with $\left|g\right|>k$ should be considered to be dubious regardless of the method of solution. Because it is difficult to extrapolate the strongly damped aerodynamic forces from the undamped aerodynamic forces via the piecewise quadratic interpolation.

Fig. 5The solutions of PQI method with mode tracking and PK method

a)$v$-$g$ plot

b)$v$-$f$ plot

Table 2 gives the real and imaginary parts of the eigenvalues ($p=g+ik$) where $\left|g\right|$ exceeds $k\text{.}$ The variation of the imaginary part of the fourth diagonal element of the generalized aerodynamic forces with the increase of $g$ is shown in Fig. 6. For the case of $g=$0, the aerodynamic forces are computed directly byusing the doublet-lattice method. For the case of $g\ne$ 0, the aerodynamic forces are computed via Eq. (2). For the large damping values such as $\left|g\right|>k$, there is a comparatively large deviation for the generalized aerodynamic forces. Therefore, for the large damping values, the aerodynamic forces should be computed directly via other methods, such as the computational fluid dynamics (CFD) method.

3.2. Flutter analysis of the advanced fighter aircraft

To verify the accuracy of the present PQI method with mode tracking in the practical engineering application, the flutter boundary of an Advanced Fighter Aircraft (AFA) with closely spaced modes is illustrated. Fig. 7 shows the configuration of the AFA model with control surfaces. It is a test bed for AeroServoelastic (ASE) analysis, and each wing of the AFA model has four control surfaces. The green-shaded areas of the AFA model denote the Leading-Edge Inboard (LEI) and Leading-Edge Outboard (LEO) control surfaces. The yellow-shaded areas denote the Trailing-Edge Inboard (TEI) and Trailing-Edge Outboard (TEO) control surfaces, respectively. The red–shaded area illustrates the tip ballast store of the AFA model. The finite element model of AFA is modeled for structural design optimization. The finite element model has an all movable horizontal tail and four control surfaces on each wing. The structural model consists of 1276 grid points and 4449 elements, and has 3761 free degrees of freedom with symmetric boundary conditions. Fig. 8 shows the finite element model of AFA, and more details on this model can be found in Ref [11-13]. In this study, the first 12 structural modes are chosen for flutter analysis. Fig. 9 shows the natural frequencies of the first 12 structural modes. As shown in the figure, the natural frequencies of the sixth and seventh modes are 23.149 Hz and 23.919 Hz, respectively. The frequency interval is merely 0.77 Hz. Moreover, the natural frequencies of eighth and ninth modes are 33.547 Hz and 33.818 Hz, respectively. The frequency interval of these two modes is only 0.271 Hz. Thus, this figure demonstrates that some of the modes are very close with each other. In this case, the closely spaced modes may induce the mode tracking problem during the flutter analysis.

Fig. 6Generalized aerodynamic forces vs. reduced frequency at different g

Fig. 7Configuration of the generic advanced fighter aircraft

Fig. 8Finite element model of AFA

Fig. 9Natural frequencies of the generic advanced fighter aircraft

For the present test case, the reduced frequency range is 0.01-0.8 and the frequency range is divided into 60 segments to construct the piecewise quadratic functions. The range of the flow speeds for flutter analysis is from 120 m/s to 360 m/s and the step size is 5 m/s. Fig. 10 shows the $v$-$g$ and $v$-$f$ diagrams by the PQI method without mode tracking. As shown in the figure, the mode tracking problem occurs between the fifth and sixth modes and this phenomenon can also be found among mode 8, mode 9, and mode 10. From this figure, one can find that, for the AFA mode with closely spaced modes, the mode tracking problem is more serious. Fig. 11 shows the computed results with mode tracking scheme. As shown in this figure, by using the present mode tracking scheme, all of the structural modes can be tracked very well.

Fig. 10The solutions of the PQI method without mode tracking

a)$v$-$g$ plot

b)$v$-$f$ plot

Fig. 11The solutions of the PQI method with mode tracking

a)$v$-$g$ plot

b)$v$-$f$ plot