2.1 Coupled dynamic model of the exciter system and test object

The structural dynamics of the test object can be expressed as

(1) \begin{equation}{{\textbf{\textit{M}}}\ddot{{\textbf{\textit{x}}}} }+ {\textbf{\textit{C}}}\ddot{{\textbf{\textit{x}}}}+{\textbf{\textit{K}}}\ddot{{\textbf{\textit{x}}}}={\textbf{\textit{f}}}\end{equation}

wherein M represents the mass matrix, C represents the damping matrix, K represents the stiffness matrix and f represents the force vector. Expanding the displacement according to the mode shape, we have

(2) \begin{equation}{\textbf{\textit{x}}}=\boldsymbol{\Phi q}\end{equation}

wherein Φ is the mode shape matrix and q is a generalised coordinate. Using equation (2) to apply a coordinate transformation, equation (1) becomes

(3) \begin{equation}\overline{{\textbf{\textit{M}}}}\ddot{{\textbf{\textit{q}}}} + \overline{{\textbf{\textit{C}}}}\ddot{{\textbf{\textit{q}}}}+ \overline{{\textbf{\textit{K}}}}{{\textbf{\textit{q}}}} = \overline{{\textbf{\textit{f}}}}\end{equation}

where $\boldsymbol{\overline{M}=\boldsymbol{\Phi}^{T}M\boldsymbol{\Phi}}$, $\boldsymbol{\overline{C}={\boldsymbol{\Phi}}^{T}C\boldsymbol{\Phi}}$, $\boldsymbol{\overline{K}={\boldsymbol{\Phi}}^{T}K\boldsymbol{\Phi}}$ and $\boldsymbol{\overline{f}={\boldsymbol{\Phi}}^{T}f}$.

The electromagnetic exciter shown in Fig. 2 is mainly composed of an exciter moving coil, supporting spring, permanent magnet, shell and base. The permanent magnet, shell and base are fixed, while the moving coil is connected to the shell through the supporting spring. Its operating principle is as follows: an alternating current is introduced into the coil to generate a varying electromagnetic field; the magnetic field interacts with the permanent magnet to drive the moving coil, causing it to vibrate up and down.

Figure 2. Diagram of an electromagnetic exciter.

According to the operating principle of the exciter, its mechanical part can be regarded as a spring–mass system with a single degree of freedom^{(Reference Zhang, Zhang, Li, Ren and Han19)}, expressed by the following mathematical model:

(4) \begin{equation}m\ddot{x}^{\prime}+c\dot{x}^{\prime}+kx^{\prime}=f^{\prime}+f^{\prime\prime}\end{equation}

where m is the mass of the exciter moving coil plus the excitation rod (here, the excitation rod is assumed to be a rigid body), c is the damping of the exciter, k is the stiffness of the exciter support, f ^′ is the reaction force of the excitation rod acting on the test piece and f ^′′ is the electromagnetic force. The electromagnetic force satisfies the formula

(5) \begin{equation}f''=BlnI=K_{f}I\end{equation}

where B is the strength of the magnetic field, l is the length of the single-turn coil, n is the number of turns in the coil and I is the input current.

The dynamic equation for the current and voltage of the exciter drive coil can be written as

(6) \begin{equation}RI+L\frac{dI}{dt}+K_{v}\dot{x}^{\prime}=V^{\prime}\end{equation}

where R is the resistance of the coil, L is the self-inductance coefficient of the coil and the voltage across the coil is V^′, $K_{v}\dot{x}^{\prime}$ while is the self-inductance voltage and K _v = K _f.

The power amplification is equivalent to the system gain in a certain frequency band, which can be modelled mathematically as

(7) \begin{equation}V^{\prime}=G_{v}V\end{equation}

where V is the input voltage of the power amplifier and G _V is the gain value.

2.2 Coupled dynamic model of the test object and exciter

In ground flutter testing, multiple exciters are required to simultaneously impose loads on the test object. Obviously, the vibration of the test object will affect the loading of the exciter force, and there is a strong coupling between the multiple exciters and the test object. For convenience of derivation, the degree of freedom of the structural dynamics of the test object is separated between excited and nonexcited points. Equation (2) can thus be written as

(8) \begin{equation}{\textbf{\textit{x}}} = \left[\begin{array}{c} x_{s}\\[5pt] x_{m}\end{array}\right] = \left[\begin{array}{c} \boldsymbol{\Phi}_{s}\\[5pt] \boldsymbol{\Phi}_{m} \end{array}\right]{\textbf{\textit{q}}}\end{equation}

where subscript s is the degree of freedom corresponding to nonexcited points, and subscript m is the degree of freedom corresponding to excited points. At the junction of the exciter and the test object, one must consider the double coordination condition for displacement and force:

(9) \begin{equation}\hspace*{14pt}{\textbf{\textit{x}}}_{m} = {\textbf{\textit{x}}}^{\prime}\end{equation}

(10) \begin{equation}{\textbf{\textit{f}}} + {\textbf{\textit{f}}}^{\prime}=0\end{equation}

Considering equations (3), (4), (6), and (7), we obtain

(11) \begin{align} \left[\begin{array}{c@{\quad}c@{\quad}c} {\overline{\textbf{\textit{M}}}}& 0&0 \\[4pt] 0 & {\textbf{\textit{m}}}&0\\[4pt] 0 & 0 & 0\end{array}\right] \left[\begin{array}{c} {\overline{\textbf{\textit{q}}}}\\[4pt] {\ddot{\textbf{\textit{x}}}}\\[4pt] \ddot{\textbf{\textit{I}}}\end{array}\right] + \left[\begin{array}{c@{\quad}c@{\quad}c} {\overline{\textbf{\textit{C}}}}& 0 & 0 \\[4pt] 0 & {\textbf{\textit{c}}} & 0\\[4pt] 0 & {\textbf{\textit{K}}}_{f} & {\textbf{\textit{L}}}\end{array}\right] \left[\begin{array}{c} \dot{\textbf{\textit{q}}}\\[4pt] \dot{\textbf{\textit{x}}^{\prime}}\\[4pt] \dot{\textbf{\textit{I}}}\end{array}\right] + \left[\begin{array}{c@{\quad}c@{\quad}c} {\overline{\textbf{\textit{K}}}} & 0 & 0 \\[4pt] 0 & {\textbf{\textit{k}}} & - {\textbf{\textit{k}}}_{\textbf{\textit{f}}}\\[4pt] 0 & 0 & {\textbf{\textit{R}}}\end{array}\right] \left[\begin{array}{c} {\textbf{\textit{q}}}\\[4pt] {{\textbf{\textit{x}}}^{\prime}}\\[4pt] {\textbf{\textit{I}}}\end{array}\right] = \left[\begin{array}{c} {\overline{\textbf{\textit{f}}}}\\[4pt] {\overline{\textbf{\textit{f}}}}^{\prime}\\[4pt] {\textbf{\textit{GV}}}\end{array}\right].\nonumber\\ \end{align}

Due to the multi-point excitation, the physical parameters associated with the exciter become diagonal matrices.

It can be shown from equations (8) and (9) that

(12) \begin{equation}{\boldsymbol{\Phi}}_{m}{\textbf{\textit{q}}} = {\textbf{\textit{x}}}'.\end{equation}

Substituting equation (12) into (11), rearranging and then considering equation (10), the following equation is obtained:

(13) \begin{align} &\left[\begin{array}{c@{\quad}c} {\overline{\textbf{\textit{M}}}}+ {\boldsymbol{\Phi}}_{m}^{T}\, {\textbf{\textit{m}}}{\boldsymbol{\Phi}}_{m} & 0 \\[6pt] 0 & 0 \\[6pt] \end{array}\right] \left[\begin{array}{c} \ddot{\textbf{\textit{q}}}\\[9pt] \ddot{\textbf{\textit{I}}} \end{array}\right] + \left[\begin{array}{c@{\quad}c} {\overline{\textbf{\textit{C}}}}+ {\boldsymbol{\Phi}}_{m}^{T}\, \textbf{\textit{c}}{\boldsymbol{\Phi}}_{m} & 0 \\[6pt] {\boldsymbol{\Phi}}_{m} {\textbf{\textit{K}}}_{f} & {\textbf{\textit{L}}} \end{array}\right] + \left[\begin{array}{c} \dot{\textbf{\textit{q}}}\\[6pt] \dot{\textbf{\textit{I}}}\end{array}\right]\nonumber\\[10pt] &\qquad + \left[\begin{array}{c@{\quad}c} {\overline{\textbf{\textit{K}}}}+ {\boldsymbol{\Phi}}_{m}^{T} {\textbf{\textit{k}}}{\boldsymbol{\Phi}}_{m} & -{\boldsymbol{\Phi}}_{m}^{T} {\boldsymbol{K}}_{f}\\[6pt] 0 & {\textbf{\textit{R}}}\\[6pt] \end{array}\right] \left[\begin{array}{c} \textbf{\textit{q}}\\[6pt] \textbf{\textit{I}}\end{array}\right] = \left[\begin{array}{c} 0\\[6pt] G{\textbf{\textit{V}}}\end{array}\right]\end{align}

For brevity, let $\boldsymbol{\tilde{M}=\overline{M}+{\boldsymbol{\Phi}}^{T}_{m}m{\boldsymbol{\Phi}}_{m}}$, $\boldsymbol{\tilde{C}=\overline{C}+{\boldsymbol{\Phi}}_{m}^{T}c{\boldsymbol{\Phi}}_{m}}$ and $\boldsymbol{\tilde{K}=\overline{K}+{\boldsymbol{\Phi}}_{m}^{T}k{\boldsymbol{\Phi}}_{m}}$. Equation (14) is then obtained from the first set of equations (13),

(14) \begin{equation}\boldsymbol{\ddot{q}}= - \boldsymbol{\tilde{M}}^{-1}\boldsymbol{\tilde{C}\ddot{q}}-\tilde{\textbf{\textit{M}}}^{-1}\boldsymbol{\tilde{K}q}+\boldsymbol{\tilde{M}}^{-1}{\boldsymbol{\Phi}}_{m}^{T}K_{f}I_{i}\end{equation}

Meanwhile, equation (13) can be further transformed into the following state-space form:

(15) \begin{equation} \left[\begin{array}{c} {\ddot{\textbf{\textit{q}}}}\\[5pt] {\dot{\textbf{\textit{q}}}}\\[5pt] {\dot{\textbf{\textit{I}}}} \end{array}\right] = \left[\begin{array}{c@{\quad}c@{\quad}c} {\tilde{-\textbf{\textit{M}}}}^{-1}\tilde{\textbf{\textit{C}}} & {\tilde{-\textbf{\textit{M}}}^{-1}\tilde{\textbf{\textit{K}}}} & {\tilde{-\textbf{\textit{M}}}^{-1}\boldsymbol{\Phi}_{m}^{T}{\textbf{\textit{K}}}_{f}}\\[5pt] \textbf{\textit{E}} & 0 &0\\[5pt] {\tilde{-\textbf{\textit{L}}}^{-1}{\textbf{\textit{K}}}_{f}\boldsymbol{\Phi}_{m}} & 0 & {-\textbf{\textit{L}}^{-1}\textbf{\textit{R}}}\end{array}\right] \left[\begin{array}{c} {\dot{\textbf{\textit{q}}}}\\[5pt] {{\textbf{\textit{q}}}}\\[5pt] {{\textbf{\textit{I}}}}\end{array}\right] + \left[\begin{array}{c} 0\\[5pt] 0\\[5pt] {\textbf{\textit{L}}^{-1}G}\end{array}\right]\textbf{\textit{V}} \end{equation}

where E represents the unit matrix. The physical quantities required for a ground flutter test system are the displacement, velocity and acceleration response of the test object, which can be solved using equation (15), then directly obtained according to equation (8). In addition, a ground flutter test system also needs to obtain the dynamic forces that the exciter applies to the test object, so the real-time MIMO controller is designed in the next section. We look again at the second set of equations in equation (11):

(16) \begin{equation}\boldsymbol{f^{\prime}}=\boldsymbol{m\Phi}_{m}\ddot{\textbf{\textit{q}}}+\boldsymbol{c\Phi}_{m}\dot{\textbf{\textit{q}}}+ \boldsymbol{k\Phi_{m}}{\textbf{\textit{q}}}-\boldsymbol{K}_{f}{\textbf{\textit{I}}}\end{equation}

Considering that the force applied by the exciter to the test object and f ^′ are reaction forces, equation (14) now yields

(17) \begin{equation} \boldsymbol{f}= \left[\begin{array}{c} \textbf{\textit{m}}\boldsymbol{\Phi}_{{m}}\tilde{\textbf{\textit{M}}}^{-1}\tilde{\textbf{\textit{C}}}-\textbf{\textit{c}}\boldsymbol{\Phi}_{m}\\[5pt] {\textbf{\textit{m}}\boldsymbol{\Phi}_{m}\tilde{\textbf{\textit{M}}}^{-1}\tilde{\textbf{\textit{K}}}-\textbf{\textit{k}}\boldsymbol{\Phi}_{m}}\\[5pt] {-\textbf{\textit{m}}}\boldsymbol{\Phi}_{m}\tilde{\textbf{\textit{M}}}^{-1}\boldsymbol{\Phi}_{m}^{T}\textbf{\textit{K}}_{f}+\textbf{\textit{K}}_{f}\\ \end{array}\right] \left[\begin{array}{c} {\dot{\textbf{\textit{q}}}}\\[5pt] {\textbf{\textit{q}}}\\[5pt] {\textbf{\textit{I}}}\end{array}\right] \end{equation}

The state-space equation consisting of equations (15) and (17) can be used directly in the design of robust MIMO dynamic force controllers.

2.3 Robust dynamic force control of the vibration exciter

Compared with numerical flutter predictions, there are more uncertainties in the coupled exciter–test model–controller–sensors coupled system. As one of the popular uncertain robust control methods, the MIMO H _∞ method can be applied to design such control laws. The mixed sensitivity method is used to design the H _∞-robust controller of the dynamic force to the excitation system, which can be transformed into a standard control problem as shown in Fig. 3.

Figure 3. Control block of H _∞-optimal control method.

The optimal control problem with the excitation system as the controlled object can be described as follows:

(18) \begin{equation} \text{min}_{K} \left||\begin{array}{c} W_{P}S\\[5pt] W_{U}KS\\[5pt] W_{G}T\end{array}\right||_{\infty} \end{equation}

where W _PS is the transfer function of the external input to the controlled output z ₁, which can be regarded as the weighted output of the error signal e, thus representing the anti-interference ability of the system. W _UKS is the transfer function of the external input to the controlled output z ₂, which can be regarded as the weighted output of the control input u, and also as the additive uncertainty output of the controlled object G. Therefore, it represents the restrictive effect on the amplitude of the control signal on the controller, as well as the robustness of the control system in the case of additive uncertainty. W _GT is the transfer function of the external input to the controlled output z ₃, which can be regarded as the multiplicative uncertainty output of the controlled object. It thus represents the robust stability of the control system in the case of multiplicative uncertainty.

The selection of the weighting functions $W_{P}(s)$, W _U(s) and W _G(s) is the key issue when designing a robust controller. The greater the amplitude of W _P(s), the higher the control accuracy of the control system. However, an excessive W _P(s) will make it more difficult to realise the control system. Generally, the amplitude is larger in the operating frequency range set by the controller, but smaller outside this range to ensure that the control system achieves good tracking characteristics. Adjusting W _U(s) will affect the control signal. The weighted function can be adjusted to prevent saturation of the control voltage. W _G(s) is the robust weighted function matrix of the system. It is determined by the non-structural uncertainty of the model, comprising the high-frequency unmodeled dynamics and the uncertainty of the model parameters. Here, note that W _P(s) and W _G(s) are generally complementary, that is, where W _P(s) is large, W _G(s) should be smaller, while W _G(s) should be larger if W _P(s) is smaller.

2.4 Fast calculation of the unsteady aerodynamic forces

For a ground flutter test system, the aerodynamic force must be calculated in real time according to the structural response of the test object. At present, in ground flutter testing systems, engineering aerodynamic models are usually applied^{(Reference Zeng, Kingsbury, Ritz, Chen, Lee and Mignolet8,Reference Xu, Wu and Yang10,Reference Song11)} , such as the dipole grid method, ZONA6, ZONA7 and so on. For an aeroelastic problem, the unsteady aerodynamic forces are calculated as

(19) \begin{equation}\boldsymbol{f}=q_{\infty}\boldsymbol{G}^{T}\boldsymbol{A}(ik)\boldsymbol{Gx}\end{equation}

where q _∞ is the dynamic pressure, G is the interpolation transformation matrix, A is the aerodynamic influence coefficient matrix, k is the reduction frequency, $k=\omega L/V$, Ω is the resonant circular frequency, L is the reference length and V is the flight speed. For convenience, the aeroelastic analysis is performed in a modal generalised coordinate system. The generalised aerodynamic model is

(20) \begin{equation} {\overline{\boldsymbol{f}}} = q_{\infty}{\boldsymbol{\Phi}^{T}\textbf{\textit{G}}^{T}}{\textbf{\textit{A}}}(ik){\textbf{\textit{G}}\boldsymbol{\Phi} \textbf{\textit{q}}}\end{equation}

For the engineering aerodynamic model, the aerodynamic influence coefficient matrix is generally expressed in the frequency domain, which must be extended to the Laplace domain and then written in state-space form^{(Reference Gao, Tan and Pu20)}. The general expression for the minimum state method rational function fitting is

(21) \begin{equation}\boldsymbol{\overline{Q}}(s)=\boldsymbol{\overline{A}}_{0}+\frac{L}{V}\boldsymbol{\overline{A}}_{1}s+\frac{L^{2}}{V^{2}}\boldsymbol{\overline{A}}_{2}s^{2}+\overline{\textbf{\textit{D}}}\left(\overline{\textbf{\textit{E}}}s-\frac{V}{L}\overline{\textbf{\textit{R}}}\right)^{-1}\overline{\boldsymbol{E}}s\end{equation}

where $\boldsymbol{Q}(s)=\boldsymbol{\Phi}^{T}{\textbf{\textit{G}}}^{T}\textbf{\textit{A}}(ik)\textbf{\textit{G}}\boldsymbol{\Phi}$, $s=ik$, $\boldsymbol{\overline{A}}_{0}$, $\boldsymbol{\overline{A}}_{1}$, $\boldsymbol{\overline{A}}_{2}$, $\boldsymbol{\overline{D}}$ and $\boldsymbol{\overline{E}}$ are matrices obtained by fitting, while $\boldsymbol{\overline{R}}$ is a diagonal matrix that is used to indicate the hysteresis state of the model caused by the hysteresis effect of the unstable airflow. The performance of the inverse Laplace transformation on equation (21) yields a time-domain calculation model whose input is the generalised response and whose output is the generalised aerodynamic force. The model is shown in Fig. 4.

Figure 4. Control block model of generalised aerodynamic forces.

4.1 Model description

The test object model is a trapezoidal aluminium plate. The length of the root chord is 2.54m, the length of the tip chord is 1.27m and the length of the semi-span is 2.54m. The support structure is located at one-third of the root of the model. The structural finite element model is shown in Fig. 5. The plate has a total number of 100 nodes and 81 shell elements. The support structure has a total number of nine nodes and nine beam elements. The support end constrains all six degrees of freedom, while the structural nodes constrain one and six degrees of freedom.

Figure 5. Finite element model of flat plate.

The parameters for the excitation system model include the mass of the exciter moving coil m, the damping c, the stiffness k, the electromagnetic constant k _f, the coil resistance R, and the inductance L as shown in Table 1. The model parameters refer to the selected exciters applied in the experimental DWT system.

4.2 Finite measurement/excitation point optimisation

The effective independent method is used to optimise the arrangement of the finite number of points for the above-described model. The measurement points are taken to be the same as the excitation points. The optimal arrangements for the cases with two to five points are shown in Fig. 6. The Fisher matrices were constructed to determine the excitation and measurement points. Because there are more than 600 Degrees of Freedom (DOFs) in this case, the Fisher matrix is very large, with dimension of more than 600 by 600. To simply the description of the designed Fisher matrix, only the largest few values of the Fisher matrix are given; For example, the matrix values of the five vertical DOFs of the five points P1 to P5 optimised in Fig. 6(d) were 0.9753, 0.9521, 0.9194, 0.8294 and 0.3515, respectively.

Figure 6. Locations of optimised points for different selected models.

After the construction of the Fisher matrix, the optimised points are then applied to the modal response estimates so that the aeroelastic analysis can be carried out. Here we used the benchmark example of NASTRAN/ZONA6 to verify the effectiveness of the proposed method for DWT. In the numerical simulation, the ZONA6 aerodynamic model was used. The frequency-domain aerodynamic model is transformed into a time-domain model by using the minimum state method. The critical flutter speeds were calculated at several Mach numbers by NASTRAN and using the proposed method with the model including five optimisation points. A comparison of the flutter speeds predicted at different Mach numbers when using both NASTRAN and the presented method is presented in Fig. 7. It seems that the flutter speed predicted by the proposed generalised force equivalence-based method agrees well with those obtained using NASTRAN. This finding indicates that the generalised force equivalence-based flutter prediction method is accurate enough to be applied as the flutter calculation solver in the DWT system.

Figure 7. Flutter speeds predicted by NASTRAN and the proposed method at different Mach numbers.

To determine the effect of the number of measurement/excitation points on the flutter speed calculation, the critical flutter speed was calculated for cases with different numbers of measurement/excitation points. The predicted errors are shown in Fig. 8, revealing that the error on the critical flutter speed obtained by the generalised force equivalent method is very small. The maximum error is less than 1% even only for the case with only two points. As the number of points is increased, the error decreases very fast and the result approaches the predicted theoretical value. From the structural dynamics perspective, the number of measurement/excitation points can theoretically be equivalent to the number of modal equations. For the aeroelastic model shown in Fig. 9, flutter is mainly caused by the first-order bending and second-order torsional modes. The first two modes can be used to describe the dominant flutter behaviour to some extent, obtaining good results even if only two points are used. As the number of points is increased, the available modal information will become more complete, so that the calculation results will approach the theoretical value. From Fig. 8, it can be concluded that, for the aeroelastic model shown in Fig. 5, the predicted errors become very small but reduce very slowly for numbers above four. The use of four points is thus a relatively good choice for DWT test systems.

Figure 8. Errors on flutter speed predicted using different numbers of sensors and exciters.

Figure 9. Aeroelastic model with measurement noise.

The discussion above considered a relatively ideal state, that is, where the modal vibration mode, measurement data, etc. all have ideal values without interference from noise. However, in real-life engineering applications, due to the influence of measurement noise and the modal identification method, certain deviations from the measurement point data and mode shape will occur. To validate the performance of the proposed method in an uncertain environment, a random mode noise interference source was introduced into the modal vibration mode data and the measurement point data in the flutter calculation procedure. The aeroelastic model with measurement noise is shown in Fig. 9. Its accuracy at Mach number 0.8 was analysed under different magnitudes of interference.

Table 2 presents the flutter speeds predicted by the generalised force equivalence-based modelling method for the aeroelastic model shown in Fig. 5 with different numbers of points and different magnitudes of noise interference. A reference flutter speed of 259 m/s is predicted by NASTRAN without noise interference at a Mach number of 0.8. As seen from Table 2, the theoretic approach using only two points can already achieve good results, although the noise has a great influence on the accuracy. As the number of points is increased, the anti-interference performance is enhanced. The maximum error is less than 1.20% for the 10% noise level when using the four-point model, and less than 1% when using the five-point model. These result indicate that the generalised force equivalence-based aeroelastic modelling method achieves good accuracy and sufficient robustness, which is very important for practical DWT systems.

Table 1 Model parameters of exciter

Table 2 Flutter speed results for different numbers of points and measurement noise levels (Reference value 259, predicted by NASTRAN without noise)

From the viewpoint of numerical simulations, more points could be selected. However, the use of more exciters requires higher-order controllers and more complex hardware systems, while the resulting increase in prediction accuracy is not so obvious. The use of an excessive number of excitation/measurement points can also cause inconvenience when performing real tests. It will lead to an increase in the order of the aerodynamic calculation model and also increase the calculation time delay of each closed loop. On the other hand, it is also necessary to increase the corresponding measuring equipment, excitation equipment and control channel, thus making the experiment more complicated. Based on comprehensive consideration of all these effects, a good choice is to use four measurement/excitation points, both to achieve good accuracy and to avoid the requirement for expensive experimental test systems for engineering applications.

4.3 Control performance of the excitation system

To analyse the influence of the accuracy of the exciter model parameters on the control system, the feedback control system model shown in Fig. 10 was constructed. In this case, the number of excitation/measurement points is four.

Figure 10. Schematic of feedback control system model for excitation system.

When designing the controller, W _P selects an eighth-order Butterworth bandpass filter with a frequency band of 5–20Hz and gain of 10. W _U takes a constant gain of 0.05, and W _G takes a constant gain of 0.1. Measurement noise interference of 10% is applied in the feedback blocks. A sinusoidal reference signal is applied to one channel, while the input reference signals from the other channels are zero. The results of the control simulation are shown in Fig. 11.

Figure 11. Output force signals of different points to input control signals of different exciters.

As seen from Fig. 11, only when a sinusoidal reference signal is applied to a certain channel can the output signal of the channel track the input reference signal well. Meanwhile, the rest of the channel outputs are close to zero, which indicates that the input of the current channel has no effect on the output of the other channels, achieving decoupling between the different channels. Moreover, the controller exhibits certain robustness, and the tracking control of the force signal and the decoupling between the different channels are achieved even in the presence of feedback noise interference. We further analyse the influence of the model parameters on the control accuracy. The system performance indicators can be calculated from the following equation based on the error signal:

(29) \begin{equation}{\textbf{\textit{error}}}=\boldsymbol{\int}^{\boldsymbol{\infty}}_{0}|\boldsymbol{e(t)}|\boldsymbol{dt}\end{equation}

In the simulation process, the variation of the system index with a certain parameter can be obtained by changing the value of a certain exciter parameter, such that sensitivity of the parameter near its nominal value is obtained; the sensitivity of the exciter parameter is obtained as shown in Fig. 12.

Figure 12. Sensitivities of exciter model parameters.

It is seen from this figure that the electrical parameters generally have a greater influence on the control performance than do the mechanical parameters such as the mass, damping and stiffness. In particular, the electromagnetic constant k _f and resistance R show the highest sensitivity, and the accuracy of these two parameters has the greatest influence on the performance of the control system. Therefore, when constructing a ground flutter test system with the proposed method for the controller design, the measurement accuracy of the electromagnetic constant k _f and resistance R should be improved to reduce the influence of parameter uncertainty on the control accuracy.

4.4 System simulation performance

Having established and validated the accuracy of each sub-system, the whole system co-simulation analysis can be carried out. The calculation process is shown in Fig. 13.

Figure 13. System-level simulation model for whole DWT system.

In this system-level simulation process, the influence of noise on the generalised modal coordinate response estimation and the time delay effect of the excitation force real-time control are considered, while the error on the critical flutter speed is also analysed (Fig. 14). The noise level indicates the measurement noise levels applied to the sensors.

Figure 14. Error in cases with different time delays and measurement noise levels.

Figure 14 confirms that the influence of the noise and the control delay of the generalised modal coordinate response estimation on the critical flutter speed cannot be ignored. Regarding the control delay, the existence of a time delay will cause the critical flutter speed to become smaller. Especially when the time delay increases to a certain extent, the error will suddenly become very large. Therefore, it is particularly important to pay attention to the influence of such time delay on the accuracy of ground flutter test systems; For example, real-time hardware devices and real-time algorithms are required for physical DWT systems.