Abstract
Establishing an accurate, fast, and low-risk flutter boundary prediction method is of great significance for flight vehicle design. In this paper, a ground flutter boundary prediction method (GFBP) based on experimental structural frequency response functions (FRFs) is proposed. A low-order multi-input multi-output (MIMO) aeroelastic system is established by combining the structural FRFs acquired from a ground test and the calculated unsteady aerodynamic FRFs in physical coordinates. The multivariable Nyquist criterion is used to predict the flutter boundary. A fixed-root aluminum plate wing is selected as the research model. A GFBP experiment is carried out for the wing’s normal state, leading-edge clump weight state, and trailing-edge clump weight state. The feasibility and accuracy of the proposed method are verified by comparison with theoretical flutter results, in which the errors of flutter speed and frequency in the test statistics are no more than 1.7%. In a simulation model established by the proposed method, Monte Carlo simulation is used to study the influence of deviations in the mode frequency and damping of the structural FRFs and deviations in the positions of excitation and measurement points in the ground test. The experiment and simulation results show that the proposed method can predict the flutter boundary accurately with accurate positions of excitation and measurement points, and it has good robustness to deviations in the mode frequency and amplitude of the structural FRFs.
1. Introduction
Flutter is a self-excited vibration of elastic structures that involves the interaction of elastic forces, inertial forces, and unsteady aerodynamic forces. This phenomenon is common in flight vehicles and may lead to serious flight accidents. Therefore, flutter boundary prediction using numerical analysis, physical tests, or semiphysical and seminumerical methods has been a continual serious concern for aviation researchers.
Numerical prediction of the flutter boundary is based on numerical models of structures and unsteady aerodynamics. The accuracy of flutter boundary prediction is determined by the accuracies of the structural models and the aerodynamic calculations. To improve the structural models, there have been many studies of semi-physical and semi-numerical methods based on the modal test. In the 1980s, Dreisbach [1] established a vibration and flutter boundary prediction method that can directly use the modal parameters (mode frequency, damping, and shape) identified by the modal test or update numerical structural models using these parameters. The concept of semi-physical and semi-numerical flutter analysis using experimental modal parameters has been extensively developed. Canfield et al. [2] used experimental modal parameters to optimize a finite element model (FEM). Zuñiga et al. [3] compare the influence of the FEM updated by experimental mode analysis and operational mode analysis on the flutter results. Pankaj et al. [4] and Chajec [5] directly incorporated experimental modal parameters into commercial software (NASTRAN or ZAERO) to predict the flutter boundary. In addition, there is a time-domain method for flutter analysis based on the structural state space model identified by MIMO modal test data. Zeng et al. [6] used this time-domain subspace identification method to establish a multi-input and multi-output (MIMO) structural state space model, which they combined with an unsteady aerodynamic state space model to perform time-domain flutter analysis. Current semi-physical and semi-numerical methods and equivalent model identification methods are used to establish the modal parameter model or state space model of the structure, which introduce identification errors in addition to the errors in the experimental data. These cumulative errors have a great impact on the accuracy of flutter boundary prediction.
For physical tests, in addition to the traditional wind tunnel test and flight test based on scaled or real models, a ground flutter simulation test (GFST) method, also called the dry wind tunnel method, based on unsteady aerodynamic real-time loading has attracted attention in recent years [7]. Distributed unsteady aerodynamic forces are calculated from structural vibration data obtained at several measurement points (MPs) as a result of excitation at several excitation points (EPs). The concentrated forces are loaded on the real structure in real time through force loading equipment. Zeng et al. [7] applied the GFST method to a rectangular flat wing and a long straight wing based on the mixed sensitivity force control. However, with high-order control, there is a restriction on the extent to which the number of EPs can be increased, making it difficult to apply this control method to complex models. Xu et al. also conducted some studies on improving the optimization method for EPs and MPs [8] and on new methods for the variable structure control [9] and the inverse dynamic model control [10] in force controller design. The orders of these force controllers are lower than control, but it is still difficult to extend them to complex models with more EPs. Hu et al. [11], Wang and Fan [12], and Zhang et al. [13] have made contributions to the GFST method in terms of condensation of aerodynamic forces, force controller design, and test error interference. For force loading, instead of an electromagnetic exciter, Yun and Ham [14, 15] used a direct drive linear actuation (DDLA) motor controlled by a simple inverse dynamic method. These studies show that the most serious problem with GFST lies in the design of the force controller. Force control errors have a significant influence on the accuracy of flutter boundary prediction. Further engineering application of this method, based as it is on a highly coupled multi-input multioutput (MIMO) system, is restricted by difficulty in realizing accurate force control.
In this study, a ground flutter boundary prediction (GFBP) method based on the structural frequency response functions (FRFs) acquired from a ground test is proposed. A low-order MIMO aeroelastic system is established by combining the experimental structural FRFs and the calculated aerodynamic FRFs in physical coordinates, and the multivariable Nyquist criterion is used to predict the flutter boundary. The structural FRFs acquired from the ground test of real structures contain more dynamical information and fewer identification errors than those obtained by identifying the modal parameters or the state space model. The proposed method adopts the same aerodynamic force condensation method as the GFST method, but it does not need real-time loading of aerodynamic forces, thereby avoiding the complex and difficult design of a MIMO force controller. Taking a fixed-root aluminum plate wing as the research model, the feasibility and accuracy of the proposed method are verified by comparison with theoretical flutter results. A simulation model established by the proposed method is used to study the effects of deviations in the mode frequency of the structural FRFs, deviations in the mode damping of the structural FRFs, and deviations in the positions of EPs and MPs in the test.
2. Methods
2.1. Principle of GFBP Method
Figure 1 shows a flow chart of the GFBP method, which involves the following steps: (1)Taking a real flight vehicle as the test model, a limited number of EPs and MPs are optimized(2)The structural FRFs are acquired from a ground test, which is performed by placing vibration sensors (laser displacement sensors, acceleration sensors, etc.) at the MPs and exciting each EP in turn(3)Establishing the aerodynamic model: using a frequency-domain unsteady aerodynamic calculation method and aerodynamic force condensation method, the aerodynamic FRFs in physical coordinates are established(4)From a combination of the experimental structural FRFs and the calculated unsteady aerodynamic FRFs, a low-order MIMO aeroelastic system is established, and the multivariable Nyquist criterion is used to predict the flutter boundary
The core of the GFBP method is to introduce the experimental structural FRFs into the frequency-domain flutter boundary prediction method. Compared with numerical models, the accuracy of structural dynamic characteristics is higher. In this paper, the research subject is limited to the linear flutter boundary prediction of flight vehicles, such as the linear flutter of wings, tails, and all-movable rudders. The nonlinear flutter phenomena such as limit cycle oscillations and chaos motions of nonlinear structures do not belong to the scope of this paper. In terms of aerodynamics, the proposed method uses mature frequency-domain unsteady aerodynamic calculation methods, needs to obtain the aerodynamic forces in physical coordinates, and uses the aerodynamic force condensation method to reduce the dimension of the aeroelastic system. The accuracy of the proposed method is limited by the accuracy of the aerodynamic method used.
2.2. Acquisition of Structural FRFs
Different from the modal test, the ground test here can also be named as frequency response test, which does not need to identify the modal parameters of the structure but only used to obtain the structural FRFs of the excitation forces at the EPs to the displacements at the MPs. From the perspective of structural and aerodynamic interaction, the structural FRFs represent the transmission characteristics from the unsteady aerodynamic forces at EPs to the displacements at MPs. The ground test can be carried out with commercial modal testing hardware and commercial software (e.g., LMS SCADAS data acquisition system and LMS Test.Lab software [16]). It should be noted that the locations of EPs and MPs cannot be given arbitrarily and need to be optimized for the accuracy of unsteady aerodynamic condensation.
Assuming that the numbers of EPs and MPs are and in the test, respectively, the dimensions of and are and , respectively, and the dimension of is . The structural subsystem is shown in Figure 2.
In the ground test, the excitation force is applied at each EP in turn, and the displacement at each MP is collected at the same time to obtain the structural FRFs. The excitation force is generally applied by an electromagnetic exciter, and the excitation force signal is collected by a force sensor. The commonly used excitation signals include burst random signal, linear sweep signal, and step-sine signal. The displacement signals can be collected directly by laser displacement sensors, or acceleration sensors can be used to collect acceleration signals that require frequency-domain integration for estimation of the structural FRFs.
In the signal processing, all kinds of FRF identification methods can be used. Commonly used methods include the H1 and H2 methods [17], which are used in the commercial modal testing software LMS Test.Lab, i.e., where and are the single-sided spectra of the auto correlation of input signal and output signal, respectively, and are all the single-sided spectra of the cross correlation of input and output signal, and is the simple harmonic oscillation frequency.
In the simulation model of the proposed method, the structural FRFs can be established from the dynamical function of the FEM: where , , and are the generalized mass matrix, the generalized damping matrix, and the generalized stiffness matrix, respectively, and and are the mode shapes of the MPs and EPs, respectively, which can be interpolated from the mode shape of the FEM.
2.3. Acquisition of Aerodynamic FRFs
In the proposed method, we need to obtain the unsteady aerodynamic FRFs in physical coordinates; the time-domain highly accurate computational fluid dynamics method based on Navier-Stokes equations can not be used directly. Therefore, the frequency-domain unsteady aerodynamic methods are used in this paper, which are fast, and with the development of theory, the accuracy generally meets the needs of engineering applications. For subsonic and supersonic speeds, the frequency-domain unsteady aerodynamic methods commonly employed in engineering can be used, such as the subsonic and supersonic doublet-lattice method (DLM) [18–20], ZONA6 [21], ZONA7 [22], and piston theory [23, 24]. The ZEUS developed by ZONA Technology, Inc. with Euler equation solver as the underlying aerodynamic force generator is a high accuracy method [25], which can be used in the transonic aerodynamic calculation.
It should be noted that for the proposed method, only the aerodynamic influence coefficient (AIC) matrix in physical coordinates needs to be calculated, and the generalized AIC matrix in generalized mode coordinates used in classical flutter analysis does not need to be obtained.
For an aeroelastic system with aerodynamic grids, the relationship between the aerodynamic forces at the pressure points and the displacements and their streamwise derivatives at the control points can be expressed as where is the atmospheric density; is the flight speed; and , , and all have dimensions . The AIC matrix in physical coordinates has dimensions . The dimensionless reduced frequency is equal to , where is the reference chord length. At a given flight speed and frequency , the AIC matrix can be obtained by interpolation.
In order to facilitate the implementation of the ground test, the number of the EPs and the MPs is generally much less than the number of aerodynamic grids. Aeroelastic system is the coupling of structural subsystem and aerodynamic subsystem, so it is necessary to ensure that their dimensions are consistent. Therefore, we adopt the aerodynamic force condensation method to obtain the FRFs from the displacement at the MPs in response to the aerodynamic forces at the EPs. The condensate aerodynamic subsystem is shown in Figure 3.
The spline interpolation methods commonly employed in general aeroelastic analysis can be used for the interpolation of displacement and force, including, among others, the infinite plate spline method (IPS) [26] and the thin plate spline method (TPS) [27].
and its streamwise derivation can be interpolated from : where and are the displacement interpolation matrices, with dimensions .
can be interpolated from the displacement at the EPs: where is the displacement interpolation matrix, with dimensions .
According to the principle of virtual work, the relationship between the aerodynamic forces at the pressure center and those at the EPs can be obtained as
On substituting Equations (4) and (6) into Equation (3), the aerodynamic forces are obtained as where the condensate AIC matrix can be written as
In order to apply the GFBP method in this paper, we take the DLM (subsonic and supersonic) as an example to calculate the unsteady aerodynamic FRFs, in which can be expressed as where is the area matrix of the aerodynamic grids and is the aerodynamic coefficient matrix, all of which have dimension . is the imaginary unit.
2.4. Optimization of EPs and MPs
Accuracy of aerodynamic force condensation is closely related to the number and locations of the EPs and the MPs. Therefore, it is necessary to optimize the locations of the EPs and the MPs. The optimization principle is that of minimizing the error between the interpolated mode shape and the original mode shape (the target mode shape) of the aerodynamic grids. This optimization principle is adopted because the smaller the mode shape error, the more accurate will be the aerodynamic force [10].
The relationship between the mode shape of the MPs and the mode shape (and the mode slope) of the aerodynamic control points is where the mode shape at the MPs can be interpolated from the experimental mode shape or assumed mode shape, is the interpolated mode shape of the aerodynamic control points, and is the mode slope.
The relationship between the mode shape of the EPs and that of the aerodynamic pressure center is where the mode shape at the EPs can be interpolated from the experimental mode shape or assumed mode shape and is the interpolated mode shape of the aerodynamic pressure center.
A genetic algorithm is used to optimize the locations of the EPs and MPs. The objective functions are expressed as where , , and are the target mode shapes (and mode slopes) (which can be interpolated from the experimental mode shape or assumed mode shape). According to practical engineering experience and the flutter characteristics of the research model, the mode weight coefficient matrix is manually set to represent the importance of each-order mode in the optimization.
It should be noted that the ground test of the proposed method is different from the GFST. The GFST needs to load unsteady aerodynamic forces on the real structure in real time. Therefore, the numbers and locations of the EPs and MPs are greatly limited by the model size, instrument size, and force control method. In this study, the test only needs to obtain the structural FRFs from the displacements at the MPs in response to the forces at the EPs, which imposes less stringent limitations on the number and locations of the EPs and MPs than the GFST.
2.5. Flutter Boundary Prediction
2.5.1. Low-Order MIMO Aeroelastic System
The structural FRFs (Figure 2) and the aerodynamic FRFs (Figure 3) are coupled to establish a low-order MIMO aeroelastic system, which is actually a feedback system, as shown in Figure 4. The inputs of the aerodynamic subsystem are the displacements at the MPs, and the outputs are the concentrated aerodynamic forces at the EPs. The FRFs of the aerodynamic subsystem are obtained by calculation. The inputs of the structural subsystem are the exciting forces at the EPs, and the outputs are the displacements at the MPs. The FRFs of the structural subsystem are obtained from a ground test of the real structure.
The flight speed can be regarded as the feedback gain in the MIMO system. When the flight speed is zero, the system is open-loop and stable. When the flight speed is less than the flutter speed, the system always remains stable. When the flight speed exceeds the flutter speed, dynamical instability of the system occurs.
2.5.2. Multivariable Nyquist Criterion in GFBP
In this study, the multivariable Nyquist criterion [28] is used to analyze the stability of the MIMO aeroelastic system and predict the flutter speed and frequency.
The open-loop FRF of the MIMO aeroelastic system at the flight speed is
The return difference matrix (RDF) is
When , the closed-loop FRF of the system is
From Equations (13)–(15), the following relationship can be established:
where is the closed-loop characteristic polynomial and is the open-loop characteristic polynomial.
Exchanging and is equivalent to the GFBP method, and so the following identity can be obtained:
In this identity, the dimensions of the matrix expressions inside the square brackets, i.e., , on the left- and right-hand sides are different: and , respectively. Therefore, it is possible to select the with the lower dimension according to the number of EPs and MPs in the test and thereby improves the calculational efficiency.
The Nyquist curve is a closed contour drawn by as transverses from to . Let denote the number of encirclements of the origin made by the curve . The number of encirclements of the origin made by is denoted by , which is equal to the number of closed-loop right-half-plane zeros . The number of encirclements of the origin made by is denoted by , which is equal to the number of open-loop right-half-plane poles . Referring to Equation (16), the following equation can be obtained:
A necessary and sufficient condition for the stability of the closed-loop system is that the closed-loop characteristic polynomial has no closed right-half-plane zeros (). Then, from Equation (18), the multivariable Nyquist criterion for the stability of the closed-loop system is that the number of clockwise encirclements relative to the origin made by plus the number of open-loop right-half-plane poles must be zero, i.e.,
In particular, for this MIMO aeroelastic system, because the structural subsystem is stable (the system is open-loop stable), has no open-loop right-half-plane poles (). The multivariable Nyquist criterion can then be represented as
In the proposed method, to avoid repetitive drawing of the Nyquist curve to judge whether it encircles the origin , this task is converted one of judging the minimum distance from to the origin. The flutter boundary prediction process is shown as a flow chart in Figure 5 and involves the following steps: (1)Take as given quantities, ranges of flight speed, and analysis frequency that cover the possible ranges of flutter speed and frequency(2)For each given flight speed, calculate the minimum distance from to the origin within the analysis frequency range(3)Within the flight speed range, obtain the - curve (as shown in Figure 6), judge the nonzero minimum value , and take the flight speed corresponding to this minimum value as the flutter speed , and take the corresponding frequency as the flutter frequency
To speed up the search for and , the flight speed range can be divided into two areas: a large-range rough search area and a small-range fine search area.
3. Research Model and Experiment
3.1. Model Parameters
In the flutter phenomenon of flight vehicles, the most common is the linear flutter of wings, tails, and all-movable rudders. In the numerical method, generally, the aerodynamic model can be simplified into flat plates with a certain geometry, and the structural model can also be simplified under the condition of ensuring the consistency of the mass and stiffness [29]. For such realistic problems, the fixed-root plate is a representative research model, which can reflect the structural and aerodynamic characteristics of wings, tails, and all-movable rudders. In the flutter boundary prediction of the fixed-root plate model using the GFBP method, the optimization method of EPs and MPs, the unsteady aerodynamic calculation and condensation method, and the experimental acquisition method of structural FRFs are consistent with the flutter boundary prediction of real flight vehicles. The research model used in this paper is shown in Figure 7, a fixed-root aluminum plate wing, with a thickness of 5 mm, a chord length of 400 mm, and a half-span length of 500 mm. The states we study include the normal state, the leading-edge clump weight state (leading state), and the trailing-edge clump weight state (trailing state). The barycenters of the leading-edge and trailing-edge clump weights are located at the coordinates (20, 480) mm and (380, 480) mm , respectively, and the weight is 95.7 g. It is a typical bending-torsional flutter model with coupling of the first bending mode and the first torsional mode. The FEM is established by grids of quadrilateral structural elements, which are used for theoretical flutter boundary prediction via method and establishing the simulation model of the proposed method. In order to improve the accuracy of the results of method, the FEM is updated by the hammering modal test values. The three lowest-order mode frequencies of the three states are listed in Table 1. It can be seen that the FEM values are in good agreement with the test values, so the theoretical flutter results can be used as the standard results of the model.
In order to verify the feasibility and accuracy of the proposed method, the same unsteady aerodynamic calculation method as that used in the theoretical method is adopted. The difference is that the condensate unsteady aerodynamic FRFs in physical coordinates need to be used in the proposed method, while the generalized AIC matrix in generalized mode coordinates is used in the theoretical method. The aerodynamic models for the three states are established by the subsonic and supersonic DLM, in which the aerodynamic grids are , the atmospheric density is taken as 1.225 kg/m3, and the Mach number is set differently. For the normal state, the calculated Mach number includes 0.2, 0.4, 0.6, and 1.5. For the leading state and trailing state, the calculated Mach number is only 0.2. Combined with the interpolation matrices (, , and ), the condensate AIC matrix in physical coordinates is computed for a selected range of reduced frequencies () and then expands by interpolation to a large number of frequencies () as required for the flutter boundary prediction.
3.2. Locations of the EPs and MPs
A genetic algorithm is used to optimize the locations of the EPs and the MPs, which are the same in the three states. Considering both the implementation efficiency of the ground test and the prediction accuracy of the flutter boundary, the numbers of EPs and the MPs are both set to four (, ). The wing model is divided into structural grids. Along the flow direction, the wing tip nodes are numbered 1–21, the root nodes are numbered 526–546, and the other nodes are numbered from small to large in this order.
Take the non-edge nodes as the search area. Considering that this is a bending-torsional flutter with coupling of the first bending mode and the first torsional mode, therefore, only the three lowest-order modes are considered, and the weight coefficient is set as . According to the optimization objective functions proposed in Equation (12), the locations of the EPs and MPs optimized by the genetic algorithm are shown in Figure 8. The numbers of EPs (E1–E4) are 390, 237, 248, and 52, and the numbers of MPs (M1–M4) are 107, 32, 76, and 125.
It should be noted that for such a fixed-root plate, a number of virtual EPs and MPs can be selected at the root of the fixed support, which can effectively improve the accuracy of the interpolation matrices , , and . This improvement makes the calculation of aerodynamic FRFs more accurate and eventually provides more accurate flutter boundary prediction results. Since the mode shape of the virtual EPs and MPs is zero, the corresponding structural FRFs are equal to zero according to Equation (12), and they do not need to be considered when establishing the MIMO aeroelastic system. When using the multivariable Nyquist criterion, it is only necessary to establish the matrix and the matrix .
Figure 9 shows a comparison of the results for the three lowest-order mode shape (and mode slope). The interpolated mode shapes and mode slopes (, , and ) are in good agreement with the target modes. The first bending mode and the first torsional mode correspond well to the target modes because these modes are relatively simple and the weight coefficient is set to 100%. For the second bending mode, because this mode is relatively complex and the weight coefficient is only set to 5%, the interpolated errors are relatively obvious, but this mode has little influence on the flutter results.
(a) and
(b) and
(c) and
3.3. Implementation of the Ground Test
The structural FRFs of the plate wing were obtained from a ground test. An electromagnetic exciter was used to excite each of the four optimized EPs, and the acceleration signals of the four MPs were collected at the same time. The FRF from the excitation force signal of a single EP to the displacement signal of a single MP was obtained by frequency-domain integration and was then combined into the structural FRFs with dimensions .
Photographs of the test implementation are shown in Figure 10. Figure 10(a) shows the four EPs excited in the normal state, and Figure 10(b) shows single-point excitation in the leading and trailing states. A TIRA S 51110-M electromagnetic exciter was used, the force signals were measured by a PCB 208C02 force sensor, and the acceleration signals were collected by PCB 333B30 acceleration sensors. All measurements and estimation of the FRFs were accomplished by the LMS SCADAS III data acquisition system and the LMS Test.Lab software package.
(a) Normal state
(b) Leading state (left) and trailing state (right)
In order to verify the quality of the experimental structural FRFs, in the normal state, the theoretical values were calculated according to Equation (2), the structural characteristics were obtained from the FEM, and the damping ratio of the three lowest-order modes was set to 0.5%.
When obtaining the structural FRFs, a burst random excitation signal was used. As Figure 11 shows, the length of the time-domain signal of one excitation was 8 s. The high-quality structural FRFs were obtained by the H1 method and averaging 20 times. Figure 12 shows a comparison of theoretical and experimental values of the structural FRFs ( and ). It can be seen that the ground test can provide high-quality structural FRFs. In the locally enlarged view at the first bending mode, there are deviations in the mode frequency and amplitude between the experimental and theoretical values, which are mainly due to the additional stiffness of the exciter and the additional mass of the sensors. The amplitudes of the FRFs are also affected by the difference in damping between the experimental and theoretical models, deviations in the positions of the EPs and MPs, and the method for estimating the FRFs. The obvious frequency difference between and shows that the additional stiffness caused by the exciter is different at different positions.
(a)
(b)
4. Results and Discussion
4.1. Prediction Results for Data from a Single Ground Test
Combined with the experimental structural FRFs and the unsteady aerodynamic FRFs, the MIMO aeroelastic system is established, and the multivariable Nyquist criterion is used to predict the flutter boundary. The GFBP results and theoretical results are shown in Table 2. It can be seen that the flutter boundary prediction results for the three states at different Mach numbers are in good agreement with the theoretical values. The experimental structural FRFs acquired from a single ground test are suitable for flutter boundary prediction at different Mach numbers.
For the normal state, the GFBP curves at Mach 0.2 are shown in Figure 13. Referring to the theoretical flutter results, the analysis speed range is 248–253 m/s, the speed step is 0.1 m/s, and the analysis frequency range is 10–60 Hz. The minimum distance between calculated from the test data and the origin is at the flutter speed 250.5 m/s. Figure 13(b) is the Nyquist curve at this critical speed.
(a) -
(b) Nyquist curve at
A further illustration of the flutter results at Mach 0.2 for the three states is shown in Figure 14. The GFBP method can accurately reflect the trends of variation of flutter speed and frequency with the clump weight located in different places. The errors in the flutter speed for the normal, leading, and trailing states are -0.44%, 1.07%, and 0.93%, respectively, while the errors in the flutter frequency are -1.42%, 0.59%, and 0.80%, respectively.
(a) Flutter speed
(b) Flutter frequency
4.2. Dispersion of Prediction Results with Data from Multiple Ground Tests
To avoid the errors and uncertainties of a single test, six groups of burst random signals with different excitation force peaks were used in the ground test. In this way, the test data from the four EPs can be combined to obtain 1296 groups of structural FRFs. Similarly, there are 1296 groups of GFBP results at Mach 0.2.
For the errors between the GFBP results of the 1296 group compared with the theoretical results, the box diagram shown in Figure 15 can be drawn, in which the maximum and minimum values of the errors are marked. It can be seen that for the three states, the results of multiple tests with different peak excitation forces fluctuate within a small range. The maximum errors in the flutter speed compared with the theoretical values for the three states are -0.60%, 1.65%, and 1.30%, respectively, while the maximum errors in the flutter frequency are -1.65%, 1.43%, and 1.03%, respectively. The GFBP results with data from multiple ground tests are stable and reliable.
(a) Flutter speed errors
(b) Flutter frequency errors
4.3. Influence of Deviations in the Structural FRFs
In the ground test, there are deviations in the mode frequency and amplitude between the experimental and theoretical FRFs. During the implementation of the test, deviations in the actual positions of EPs and MPs will also lead to errors in the FRFs. In the simulation system for the normal-state research model established by the proposed method, the mode frequency, the mode damping (reflecting amplitude deviation), and the positions of EPs and MPs are each randomly perturbed, and Monte Carlo simulation is used to count the flutter boundary prediction results. The unbiased GFBP results for the normal-state simulation model at are shown in Figure 16, namely, the - curve representing the minimum distance from to the origin and the Nyquist curve at the flutter speed. For the unbiased results, the flutter speed is 250.3 m/s, and the flutter frequency is 30.90 Hz.
(a) -
(b) Nyquist curve at
4.3.1. Deviations in the Mode Frequency
To get close to the real situation of the ground test, considering that the mode frequencies obtained at different EPs are affected differently by the additional stiffness of the exciter and the additional mass of the sensors, the frequency of each-order mode obtained at each EP is subjected to independent random perturbations. The random perturbation of each-order mode frequency is given by where is the nominal mode frequency of the normal state, is the perturbation mode frequency, is a random variable satisfying a Gaussian distribution, and represents the mode frequency perturbation range of the 95% confidence interval. Two perturbation ranges of 1% and 2% are considered in the analysis; here, 2% is the maximum range of mode frequency deviation in the test of the plate wing.
Using Monte Carlo simulation, 4000 groups of samples were randomly generated, and the GFBP results of each group were counted by the kernel smoothing method. The probability density function (PDF) of the flutter speed is shown in Figure 17. When the frequency deviation conforms to a Gaussian distribution, the predicted flutter speed also conforms to a Gaussian distribution. For the 1% perturbation range, the mean flutter speed is 250.3 m/s, the 95% confidence interval is 247.6 to 253.0 m/s, and the error in the 95% confidence interval relative to the unbiased state is . For the 2% perturbation range, the mean flutter speed is 250.3 m/s, the 95% confidence interval is 244.8 to 255.8 m/s, and the error in the 95% confidence interval relative to the unbiased state is .
The analysis of the frequency deviation shows that the maximum possible frequency deviation range in the test has little influence on the GFBP results. The proposed method has good robustness and can give accurate flutter boundary prediction results.
4.3.2. Deviations in the Mode Damping
The deviation of mode damping is mainly used to reflect the amplitude deviation of the structural FRF. Considering that the amplitude deviation at different EPs is different, to be closer to the real situation, the damping ratios of the three lowest-order modes obtained at each EP are subjected to independent random perturbations. The nominal damping ratio is 0.5%, and the perturbation damping ratio is given by where is a random variable satisfying a Gaussian distribution and represents the mode damping perturbation range of the 95% confidence interval. Two perturbation ranges of 20% and 60% are considered in the analysis.
Using Monte Carlo simulation, 4000 groups of samples were randomly generated, and the GFBP results of each group were counted by the kernel smoothing method. The PDF of the flutter speed is shown in Figure 18. It can be seen that the smaller the perturbation of the damping ratio, the smaller is the variance of the flutter speed. For the 20% perturbation range, the mean flutter speed is 250.3 m/s, and the 95% confidence interval is 249.8 to 250.8 m/s. For the 60% perturbation range, the mean flutter speed is 250.3 m/s, the 95% confidence interval is 248.7 to 251.9 m/s, and the error in the 95% confidence interval relative to the unbiased state is . This result indicates that for this kind of explosive flutter, damping has little effect on the flutter speed. When the GFBP method is used to predict the flutter boundary, there is a large allowable deviation in the damping or, more directly, the amplitude of the structural FRFs.
4.3.3. Deviations in the Positions of the EPs and MPs
In the test, owing to artificial measurement errors, the actual positions of the EPs and MPs may deviate from the theoretical values. To study the influence of these deviations on the GFBP results, the position coordinates of the EPs and MPs are randomly perturbed. The perturbation range is a circular area, with center at the theoretical coordinates and radius 5 mm (1.0% of the maximum characteristic length of the model). The perturbed coordinates are used to calculate the structural FRFs, while the aerodynamic FRFs are still calculated using the theoretical coordinates.
Four thousand groups of GFBP results were randomly generated and counted by the kernel smoothing method. The PDF of the flutter speed is shown in Figure 19 and can be seen to conform to a Gaussian distribution. The mean value of the flutter speed is 251.2 m/s, the 95% confidence interval is 232.5 to 269.9 m/s, and the error relative to the unbiased state is -7.1% to 7.8%. It can be seen that the position deviation has a great influence on the flutter boundary prediction results. The reason is that in the calculation of the aerodynamic FRFs, the theoretical coordinates of EPs and MPs are used to construct the interpolation matrices , , and , while the obtained structural FRFs deviate from the unbiased state.
Considering the different influences of the EPs and MPs in the aerodynamic calculation, their positional coordinates are perturbed separately, although the perturbation range remains unchanged, and 4000 groups of samples are still generated randomly. The PDFs of the flutter speed in the two perturbation cases are shown in Figure 20. The mean flutter speed of the EP position deviation is 250.4 m/s, the 95% confidence interval is 245.1 to 255.7 m/s, and the error relative to the unbiased state is -2.1% to 2.2%. The mean flutter speed of the MP position deviation is 251.2 m/s, the 95% confidence interval is 233.2 to 269.1 m/s, and the error relative to the unbiased state is -6.8% to 7.5%, which is close to the simultaneous perturbation results for the EPs and MPs. This result shows that compared with the positional accuracy of the EPs, that of the MPs has a greater impact on the GFBP results.
(a) Deviation in the positions of Eps
(b) Deviation in the positions of Mps
The study shows that when obtaining the structural FRFs from the test, attention should be paid to the positional accuracy of EPs and MPs, especially that of the MPs. In a ground test of a real wing, the accuracy of the positional coordinates can be guaranteed by the high-precision dimensional measuring equipment.
5. Conclusions
In this paper, a ground flutter boundary prediction method based on the structural FRFs acquired from a ground test has been proposed. Combined with experimental structural FRFs and calculated unsteady aerodynamic FRFs in physical coordinates, a low-order MIMO aeroelastic system has been established, and the multivariable Nyquist criterion has been used to predict the flutter boundary. The ground test can be carried out with commercial modal testing hardware and commercial software, and the accuracy of frequency-domain unsteady aerodynamic calculation methods commonly used in engineering generally meets the needs of industrial applications, which make this method very attractive for engineering applications. Compared with the previous flutter boundary prediction methods, the proposed method has several contributions. Firstly, compared with the numerical method base on the numerical structural model, the structural FRFs obtained from the ground test of the real structure can more accurately reflect the real dynamic characteristics. Secondly, the experimental structural FRFs have fewer identification errors than the modal parameters and state space model identified from test data. Finally, compared with the GFST method, there is no risk of instability and no need to use the complex and difficult MIMO force controller in the ground test, which makes the proposed method safe and easy to implement, and avoids the influence of force control error on the flutter boundary prediction results.
A GFBP experiment of an aluminum plate wing has been carried out. The structural FRFs are acquired from ground tests with different excitation peaks. The unsteady aerodynamic FRFs are calculated by subsonic and supersonic DLM in physical coordinates. When the numbers of EPs and MPs are both four, the maximum errors of the flutter speed and frequency in the test statistics for the normal state, leading state, and trailing state compared with the theoretical values are no more than 1.7%, which verifies the feasibility and accuracy of the proposed method. The structural FRFs acquired from a ground test can be used to flutter boundary prediction at different Mach numbers.
For the common structural FRF errors in the ground test, including deviations in the mode frequency, deviations in the mode damping, and the position deviations of EPs and MPs, a simulation model of the proposed method has been established in this paper. The results of a perturbation analysis of the simulation system show that deviations in the mode frequency of the structural FRFs have limited impact on the accuracy of the proposed method. Deviations in the damping (represent the amplitude of the structural FRFs) also have little impact for this kind of explosive flutter. Deviations in the positions of MPs have a much greater impact on the flutter boundary prediction results than deviations in the positions of EPs, and so the positional accuracy of the MPs needs to be guaranteed in the ground test.
Data Availability
The numerical and experimental data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
The support of the Aeronautical Science Foundation of China (Grant No. 20200001051001) is greatly appreciated.