2.1 Identification of hub dynamic loads by vibration measurement on fuselage

For a multi-degree-of-freedom structural system, under the exciting load ${F_j}\left( \omega \right)$ at point j of the structure, the response generated at point i of the structure is ${U_i}\left( \omega \right)$ , and there is the following formula:

(1) \begin{equation}{H_{ij}}\left( \omega \right) = \frac{{{U_i}\left( \omega \right)}}{{{F_j}\left( \omega \right)}}\end{equation}

where, ${H_{ij}}\left( \omega \right)$ is the FRF between points i and j, its physical meaning is the response generated at the point i when a unit exciting load with frequency ω is applied on point j of the structure. Due to the reciprocity of linear systems and the fuselage is a linear structure, there is the following equation:

(2) \begin{equation}{H_{ij}}\left( \omega \right) = {H_{ji}}\left( \omega \right)\end{equation}

According to the FRF definition by Equation (1), when the structure is only excited by a load ${F_j}\left( \omega \right)$ , there is the following equation:

(3) \begin{equation}{U_i}\left( \omega \right) = {H_{ij}}\left( \omega \right){F_j}\left( \omega \right)\end{equation}

If the structure is excited by N loads $\left\{ F \right\} = \left\{ {{F_1},{F_2}, \cdots {F_N}} \right\}$ , according to the principle of linear superposition, there is the following equation:

(4) \begin{equation}{U_i} = {H_{i1}}{F_1} + {H_{i2}}{F_2} + \cdots + {H_{in}}{F_n} = \left[ {{H_{i1}},{H_{i2}}, \cdots ,{H_{in}}} \right]\left\{ {\begin{array}{*{20}{c}}{{F_1}}\\[3pt]{{F_2}}\\[3pt] \vdots \\[3pt]{{F_n}}\end{array}} \right\}\end{equation}

Then there is the following equation:

(5) \begin{equation}\left\{ U \right\} = \left\{ {\begin{array}{*{20}{c}}{{U_1}}\\[3pt]{{U_2}}\\[3pt] \vdots \\[3pt]{{U_n}}\end{array}} \right\} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{H_{11}}} & {{H_{12}}} & \cdots & {{H_{1n}}}\\[3pt]{{H_{21}}} & {{H_{22}}} & \cdots & {{H_{2n}}}\\[3pt] \vdots & \vdots & & \vdots \\[3pt]{{H_{n1}}} & {{H_{n2}}} & \cdots & {{H_{nn}}}\end{array}} \right]\left\{ {\begin{array}{*{20}{c}}{{F_1}}\\[3pt]{{F_2}}\\[3pt] \vdots \\[3pt]{{F_n}}\end{array}} \right\} = \left[ H \right]\left\{ F \right\}\end{equation}

where the FRF matrix $\left[ H \right]$ is a symmetric matrix.

For the deterministic vibration responses of the structure, it is assumed that the number of exciting loads to be identified is ${n_f}$ and the number of vibration measuring points on the structure is ${n_u}$ . It is required that the number of measuring points is greater than or equal to the number of exciting loads, i.e., ${n_u} \geqslant {n_f}$ ; otherwise, the solution cannot be solved. In this case, the column vector of exciting loads ${\left[ {F\left( \omega \right)} \right]_{{n_f} \times 1}}$ and the column vector of measured responses ${\left[ {U\left( \omega \right)} \right]_{{n_u} \times 1}}$ in frequency domain of the vibration system satisfy the following relationship:

(6) \begin{equation}\left[ {U\left( \omega \right)} \right] = \left[ {H\left( \omega \right)} \right]\left[ {F\left( \omega \right)} \right]\end{equation}

where the FRF matrix ${\left[ {H\left( \omega \right)} \right]_{{n_u} \times {n_f}}}$ is a square matrix when ${n_u} = {n_f}$ , and the exciting loads can be calculated by the following formula:

(7) \begin{equation}\left[ {F\left( \omega \right)} \right] = {\left[ {H\left( \omega \right)} \right]^{ - 1}}\left[ {U\left( \omega \right)} \right]\end{equation}

Equation (7) contains the inverse of the FRF matrix ${\left[ {H\left( \omega \right)} \right]^{ - 1}}$ . The ${\left[ {F\left( \omega \right)} \right]_{{n_f} \times 1}}$ obtained by Equation (7) is the column vector of exciting loads in frequency domain. To get the time history of exciting loads, it needs to take the inverse Fourier transform of the exciting loads.

In order to improve the identification accuracy and stability of the exciting loads, the number of vibration measuring points on the structure is generally more than the number of exciting loads to be identified, i.e., ${n_u} \gt {n_f}$ , then the FRF matrix is of order ${n_u} \times {n_f}$ , and the exciting loads can be calculated by the following formula:

(8) \begin{equation}\left[ {F\left( \omega \right)} \right] = {\left[ {H\left( \omega \right)} \right]^ + }\left[ {U\left( \omega \right)} \right] = {\left[ {{{\left[ {H\left( \omega \right)} \right]}^H}\left[ {H\left( \omega \right)} \right]} \right]^{ - 1}}\left[ {{{\left[ {H\left( \omega \right)} \right]}^H}\left[ {U\left( \omega \right)} \right]} \right]\end{equation}

where superscript + represents the generalised inverse of matrix and superscript H represents the complex conjugate transpose of matrix.

Hence, according to Equation (8), the exciting loads $\left[ {F\left( \omega \right)} \right]$ can be identified as long as the vibration responses $\left[ {U\left( \omega \right)} \right]$ of the structure under the exciting loads $\left[ {F\left( \omega \right)} \right]$ are measured and the FRF matrix $\left[ {H\left( \omega \right)} \right]$ are obtained by the measurements between the exciting and measuring points on the structure. For the identification of the hub dynamic loads, the exciting points of the hub dynamic loads are determined. Hence, the measured FRF matrix on the fuselage depends on the layout of the vibration measuring points (including location and direction) on the fuselage. However, the layout of the measuring points on the fuselage is uncertain generally and may be unreasonable, which may lead to the ill-conditioned FRF matrix and then the great errors of the identification of the hub dynamic loads. The unreasonable layout means that the signal-to-noise ratios on the unreasonable measuring points are very small, and the signals measured on these measuring points contain lots of noise. The lots of error information or insufficient modal information will cause the ill-conditioned FRF matrix to produce inaccurate identification of the hub dynamic loads.

2.2 Condition number of measured FRF matrix

In order to avoid the ill-conditioned or even singular FRF matrix measured on the fuselage due to unreasonable layout of the measuring points on the fuselage, a method for reasonable layout of the measuring points on the fuselage based on the condition number of the FRF matrix measured on the fuselage is proposed in this section. The condition number of a matrix is a measure of the accuracy of numerical solution of a linear system of equations. If the condition number is not too much larger than one, the matrix is well-conditioned,which means its inverse can be computed with good accuracy. If the condition number is very large, then the matrix is said to be ill-conditioned. Such a matrix is almost singular, and the computation of its inverse, or solution of a linear system of equations is prone to large numerical errors. A matrix that is not invertible has condition number equal to infinity.

For the FRF matrix $\left[ {H\left( \omega \right)} \right]$ in Equation (8), the condition number $cond(H)$ of the matrix $\left[ {H\left( \omega \right)} \right]$ can be defined as the product of the matrix norm and its inverse matrix norm, i.e.:

(9) \begin{equation}cond(H) = \left\| H \right\| \cdot \left\| {{H^{ - 1}}} \right\|\end{equation}

where the condition number $cond(H)$ of any matrix is always greater than or equal to one. The condition number of a singular matrix is infinite, and the condition number of an orthogonal matrix is one. The greater the condition number of a matrix, the more ill conditioned the matrix.

The condition number of a matrix is related to the norm selected. In this paper, the two norm of matrix is selected as the norm of the FRF matrix $\left[ {H\left( \omega \right)} \right]$ . Hence, we have the following formula for calculating the condition number $cond(H)$ of the FRF matrix $\left[ {H\left( \omega \right)} \right]$ :

(10) \begin{equation}cond(H) = \frac{{{\sigma _{\max }}(H)}}{{{\sigma _{\min }}(H)}}\end{equation}

Where ${\sigma _{\max }}(H)$ and ${\sigma _{\min }}(H)$ are the maximum and minimum singular values of the FRF matrix $\left[ {H\left( \omega \right)} \right]$ , respectively.

Figure 1. Photo of frame structure.

In practice, the unit hub dynamic loads are applied at the hub centre of the fuselage, then the vibration responses at the measuring points on the fuselage are measured, and then the FRF matrix $\left[ {H\left( \omega \right)} \right]$ between the exciting points of the hub dynamic loads and the measuring points of vibration responses are obtained according to the definition of the FRF by Equation (1). After obtaining the FRF matrix $\left[ {H\left( \omega \right)} \right]$ , the condition number $cond(H)$ of the FRF matrix $\left[ {H\left( \omega \right)} \right]$ are calculated by Equation (10), then we can judge if the FRF matrix $\left[ {H\left( \omega \right)} \right]$ is ill-conditioned according to the value of the condition number $cond(H)$ . If the value of the condition number $cond(H)$ is not much more than one, then the FRF matrix $\left[ {H\left( \omega \right)} \right]$ is not ill-conditioned, and then the layout of the measuring points corresponding to the FRF matrix $\left[ {H\left( \omega \right)} \right]$ is reasonable. After obtaining the reasonable layout of the measuring points on the fuselage, we can identify the hub dynamic loads by measuring the vibration responses at the reasonable layout of the measuring points under the excitation of the hub dynamic loads according to Equation (8).

3.1 Fuselage model structure

The fuselage model structure used for the identification of the hub dynamic loads is a dynamically similar frame structure of the fuselage cockpit floor of Z-11 helicopter. A photo of the frame structure which was vertically installed to the foundation is shown in Fig. 1. The frame structure was assembled by rectangular cross section beams and fixed to the foundation by 10 cantilever beams that were used to simulate the connection between the fuselage cockpit floor and the middle fuselage, which is the foundation in Fig.1. The material of beam is steel and its parameters are listed in Table 1. The schematic diagram of the frame structure with the loading point, boundary nodes connecting to the foundation and verifying points 1 and 2 for verifying the response calculation method by experiment is shown in Fig. 2.

Table 1. Parameters of frame structure

Table 2. Natural frequencies of frame structure

Figure 2. Schematic diagram of frame structure.

In this paper, the finite element model of the frame structure was established by using three-dimensional beam element. The simulation and measurement results of the natural frequencies of the frame structure are listed in Table 2. It can be seen from Table 2 that the errors between the natural frequencies obtained by the finite element simulation and measurement of the frame structure are all less than 3%, indicating that the frequencies between the finite element model and the actual frame structure are very close.

3.2 Identification simulation of hub dynamic loads

Helicopter rotor hub dynamic loads are three forces and three moments in three directions which are located at the hub centre and produce vibration of fuselage structure. According to Section 2, the hub dynamic loads can be identified by obtaining the FRF matrix at the reasonable layout of the measuring points on the fuselage. In this paper, the time domain vibration response of the finite element model of the frame structure under the excitation of dynamic loads is solved by using the Newmark-Beta method.

For the vibration frame structure of the fuselage cockpit floor shown in Fig. 1 with mass matrix $\left[ M \right]$ , damping coefficient matrix $[C]$ and stiffness matrix $\left[ K \right]$ , the equation of motion under the dynamic load vector $\left\{ F \right\}$ at time point $i$ is

(11) \begin{equation}\left[ M \right]\left\{ {{{\ddot d}_i}} \right\} + [C]\left\{ {{{\dot d}_i}} \right\} + \left[ K \right]\left\{ {{d_i}} \right\} = \left\{ {{F_i}} \right\}\end{equation}

The Newmark-Beta method approximately solves the velocity and displacement of a vibrating system after a time increment $\Delta t$ using the following two equations:

(12) \begin{equation}\left\{ {{{\dot d}_{i + 1}}} \right\} = \left\{ {{{\dot d}_i}} \right\} + \Delta t\left[ {\left( {1 - \gamma } \right)\left\{ {{{\ddot d}_i}} \right\} + \gamma \left\{ {{{\ddot d}_{i + 1}}} \right\}} \right]\end{equation}

(13) \begin{equation}\left\{ {{d_{i + 1}}} \right\} = \left\{ {{d_i}} \right\} + \Delta t\left\{ {{{\dot d}_i}} \right\} + {\left( {\Delta t} \right)^2}\left[ {\left( {\frac{1}{2} - \beta } \right)\left\{ {{{\ddot d}_i}} \right\} + \beta \left\{ {{{\ddot d}_{i + 1}}} \right\}} \right]\end{equation}

where the parameters β and γ are freely selectable. Generally, the value of parameter β is between 0 and 0.25, and the value of parameter γ is 0.5. For the combination of β=0.25 and γ=0.5 values, namely the averaged acceleration method, the numerical analysis is stable and the calculation results are convergent regardless of the time step size.

By Equation (13), the following equation is obtained:

(14) \begin{equation}\left\{ {{{\ddot d}_{i + 1}}} \right\} = \frac{1}{{\beta {{\left( {\Delta t} \right)}^2}}}\left( {\left\{ {{d_{i + 1}}} \right\} - \left\{ {{d_i}} \right\}} \right) - \frac{1}{{\beta \Delta t}}\left\{ {{{\dot d}_i}} \right\} - \left( {\frac{1}{{2\beta }} - 1} \right)\left\{ {{{\ddot d}_i}} \right\}\end{equation}

Substituting Equation (14) into Equation (12) yields:

(15) \begin{equation}\left\{ {{{\dot d}_{i + 1}}} \right\} = \frac{\gamma }{{\beta {{\left( {\Delta t} \right)}^2}}}\left( {\left\{ {{d_{i + 1}}} \right\} - \left\{ {{d_i}} \right\}} \right) + \left( {1 - \frac{\gamma }{\beta }} \right)\left\{ {{{\dot d}_i}} \right\} + \left( {1 - \frac{\gamma }{{2\beta }}} \right)\Delta t\left\{ {{{\ddot d}_i}} \right\}\end{equation}

At time point i+1, the equation of motion of the vibration system is

(16) \begin{equation}\left[ M \right]\left\{ {{{\ddot d}_{i + 1}}} \right\} + [C]\left\{ {{{\dot d}_{i + 1}}} \right\} + \left[ K \right]\left\{ {{d_{i + 1}}} \right\} = \left\{ {{F_{i + 1}}} \right\}\end{equation}

Substituting Equations (14) and (15) into Equation (16) yields:

(17) \begin{equation}\left[ {\tilde K} \right]\left\{ {{d_{i + 1}}} \right\} = \left\{ {{{\tilde F}_{i + 1}}} \right\}\end{equation}

(18) \begin{equation}\left[ {\tilde K} \right] = \frac{1}{{\beta {{\left( {\Delta t} \right)}^2}}}\left[ M \right] + \frac{\gamma }{{\beta \Delta t}}[C] + \left[ K \right]\end{equation}

(19) \begin{align}\left\{ {{{\tilde F}_{i + 1}}} \right\} & = \left\{ {{F_{i + 1}}} \right\} + \left[ M \right]\left[ {\frac{1}{{\beta {{\left( {\Delta t} \right)}^2}}}\left\{ {{d_i}} \right\} + \frac{1}{{\beta \Delta t}}\left\{ {{{\dot d}_i}} \right\} + \left( {\frac{1}{{2\beta }} - 1} \right)\left\{ {{{\ddot d}_i}} \right\}} \right]\nonumber\\[5pt]& \quad + [C]\left[ {\frac{\gamma }{{\beta \Delta t}}\left\{ {{d_i}} \right\} + \left( {\frac{\gamma }{\beta } - 1} \right)\left\{ {{{\dot d}_i}} \right\} + \left( {\frac{\gamma }{{2\beta }} - 1} \right)\Delta t\left\{ {{{\ddot d}_i}} \right\}} \right]\end{align}

Hence, according to Equations (17), (15) and (14), the displacement, velocity and acceleration at time point i+1 after a time increment $\Delta t$ can be solved by using the parameters β, γ, $\Delta t$ and the values of displacement, velocity, acceleration and external force at time point $i$ .

For verifying the response calculation using the Newmark-Beta method by experiment, the simulation and experiment at the verifying points 1 and 2 as shown in Fig. 2 were carried out. In the simulation, a dynamic force ${F_z}$ with frequency 25Hz and amplitude 20N and a dynamic moment ${M_y}$ with frequency 25Hz and amplitude 14Nm were applied at the loading point, respectively, and the acceleration responses at the verifying points 1 and 2 were calculated by using the Newmark-Beta method under the dynamic force and torque, respectively. In the experiment, a dynamic HEV-50 shaker was used to produce the same dynamic force as in simulation and apply to the loading point, two dynamic HEV-50 shakers with a distance of 14cm and opposite phases were used to produce the same dynamic torque as in simulation and apply to the loading point, and the acceleration responses at the verifying points 1 and 2 were measured. The simulation and experimental acceleration amplitudes at two verifying points are listed in Table 3. It can be seen from Table 3 that the largest error of the simulation with experiment is 6.22%, verifying the calculation results using the Newmark-Beta method.

Table 3. Simulation and experimental acceleration amplitudes at two verifying points

For the frame structure shown in Fig. 2, the loading point of the hub dynamic loads and six measuring points are selected as shown in Fig. 3 in which the black dot is the exciting point of the hub dynamic loads, the square black dots are the measuring points, the arrow is the measurement direction of vibration acceleration, and x, y and z refer to the lateral, longitudinal and vertical directions of the fuselage, respectively. The rotor hub dynamic loads include the three forces ${F_x},{\rm{ }}{F_y},{\rm{ }}{F_z}$ and three moments ${M_x},{\rm{ }}{M_y},{\rm{ }}{M_z}$ . However, in the simulation, the lateral force ${F_x}$ , vertical force ${F_z}$ , longitudinal moment ${M_y}$ and vertical moment ${M_z}$ were used as the hub dynamic loads. The longitudinal force ${F_y}$ and lateral moment ${M_x}$ were not used in the simulation, because the end of the frame structure is fixed on the foundation, the vibration of the frame structure under ${F_y}$ and ${M_x}$ is too small to effectively identify the dynamic loads. Hence, four dynamic loads, one loading point and six measuring points were adopted in the identification simulation of the hub dynamic loads.

Figure 3. Loading point of hub dynamic loads and measuring points on frame structure.

The unit dynamic loads were applied to the loading point and the acceleration responses on the six measuring points were computed by using the Newmark-Beta method to obtain three FRF matrices corresponding to 18, 36 and 54Hz which are the 1st, 2nd and 3rd blade passage frequencies of Z-11 helicopter, respectively. The three $6 \times 4$ order FRF matrices are denoted as $H\left( {18} \right)$ , $H\left( {36} \right)$ and $H\left( {54} \right)$ , respectively. In order to simulate the background noise, Gaussian white noise with an amplitude range of $ \pm {10^{ - 2}}m/{s^2}$ was added to the vibration response signal in time domain. According to the simulation results, the acceleration responses at measuring point 3 under the excitation of ${F_x}$ , at measuring points 1, 4 and 6 under the excitation of ${F_z}$ , at measuring point 3 under the excitation of ${M_y}$ and at measuring point 3 under the excitation of ${M_z}$ were all white noise responses. The simulated accelerations at six measuring points under unit dynamic force ${F_x}$ are shown in Fig. 4 in which the acceleration signal at measuring point 3 is actually the noise. At other measuring points and under the excitation of other dynamic loads, the signal-to-noise ratios of vibration acceleration responses were very high, which indicated that the vibration acceleration responses were large.

Figure 4. Simulated accelerations at six measuring points under unit dynamic force ${F_x}$ .

According to the applied unit dynamic loads and the vibration responses at the measuring points, the FRF matrices $H\left( {18} \right)$ , $H\left( {36} \right)$ and $H\left( {54} \right)$ of the frame structure can be obtained as follows:

(20) \begin{equation}H\left( {18} \right) = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{\rm{ - 0}}{\rm{.0058 + 0}}{\rm{.0001i}}} & {\rm{0}} & {{\rm{0}}{\rm{.024 + 0}}{\rm{.0003i}}} & {{\rm{0}}{\rm{.036 - 0}}{\rm{.0005i}}}\\[3pt]{{\rm{0}}{\rm{.0064 - 0}}{\rm{.0001i}}} & {{\rm{ - 0}}{\rm{.012}} + {\rm{0}}{\rm{.0072i}}} & {{\rm{ - 0}}{\rm{.085 + 0}}{\rm{.0025i}}} & {{\rm{ - 0}}{\rm{.040 + 0}}{\rm{.0010i}}}\\[3pt]{\rm{0}} & {{\rm{ - 0}}{\rm{.13 + 0}}{\rm{.019i}}} & {\rm{0}} & {\rm{0}}\\[3pt]{{\rm{ - 0}}{\rm{.0091}} + {\rm{0}}{\rm{.0003i}}}& {\rm{0}} & {{\rm{ - 0}}{\rm{.044 + 0}}{\rm{.0011i}}} & {{\rm{0}}{\rm{.091 - 0}}{\rm{.0026i}}}\\[3pt]{{\rm{0}}{\rm{.0052 - 0}}{\rm{.0001i}}}& {{\rm{ - 0}}{\rm{.17 + 0}}{\rm{.018i}}} & {{\rm{ - 0}}{\rm{.069 + 0}}{\rm{.0012i}}}& {{\rm{ - 0}}{\rm{.031 + 0}}{\rm{.0006i}}}\\[3pt]{{\rm{ - 0}}{\rm{.018 + 0}}{\rm{.0006i}}} & {\rm{0}} & {{\rm{0}}{\rm{.048 - 0}}{\rm{.0013i}}}& {{\rm{0}}{\rm{.16 - 0}}{\rm{.0050i}}}\end{array}} \right]\end{equation}

(21) \begin{equation}H\left( {36} \right) = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{\rm{ - 0}}{\rm{.0096}} + {\rm{0}}{\rm{.0002i}}} & {\rm{0}} & {{\rm{0}}{\rm{.026}} - {\rm{0}}{\rm{.0013i}}} & {{\rm{0}}{\rm{.036 - 0}}{\rm{.0019i}}}\\[3pt]{{\rm{ - 0}}{\rm{.015 - 0}}{\rm{.0006i}}} & {{\rm{1}}{\rm{.1}} + {\rm{0}}{\rm{.18i}}} &{{\rm{0}}{\rm{.20}} + {\rm{0}}{\rm{.0081i}}} &{{\rm{0}}{\rm{.094}} + {\rm{0}}{\rm{.0038i}}}\\[3pt]{\rm{0}} & {{\rm{0}}{\rm{.50}} - {\rm{0}}{\rm{.046i}}} & {\rm{0}} & {\rm{0}}\\[3pt]{{\rm{0}}{\rm{.019}} + {\rm{0}}{\rm{.0005i}}} &{\rm{0}} & {{\rm{0}}{\rm{.13}} + {\rm{0}}{\rm{.0041i}}} & {{\rm{ - 0}}{\rm{.20}} - {\rm{0}}{\rm{.0057i}}}\\[3pt]{{\rm{ - 0}}{\rm{.0011 - 0}}{\rm{.0004i}}} & {{\rm{ - 1}}{\rm{.1}} + {\rm{0}}{\rm{.19i}}} & {{\rm{0}}{\rm{.014}} + {\rm{0}}{\rm{.0054i}}} & {{\rm{0}}{\rm{.0084}} + {\rm{0}}{\rm{.0026i}}}\\[3pt]{{\rm{0}}{\rm{.071}} + {\rm{0}}{\rm{.0014i}}} & {\rm{0}} &{{\rm{ - 0}}{\rm{.19}} - {\rm{0}}{\rm{.0046i}}} &{{\rm{ - 0}}{\rm{.59 - 0}}{\rm{.012i}}}\end{array}} \right]\end{equation}

(22) \begin{equation}H\left( {54} \right) = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{\rm{ - 0}}{\rm{.054}} + {\rm{0}}{\rm{.0015i}}} & {\rm{0}} & {{\rm{0}}{\rm{.21}} - {\rm{0}}{\rm{.0079i}}} & {{\rm{0}}{\rm{.33 - 0}}{\rm{.011i}}}\\[5pt]{{\rm{0}}{\rm{.016}} - {\rm{0}}{\rm{.0018i}}} & {{\rm{ - 0}}{\rm{.32}} + {\rm{0}}{\rm{.021i}}} &{{\rm{ - 0}}{\rm{.22}} + {\rm{0}}{\rm{.024i}}} &{{\rm{ - 0}}{\rm{.099}} + {\rm{0}}{\rm{.011i}}}\\[5pt]{\rm{0}} & {{\rm{ - 0}}{\rm{.46}} - {\rm{0}}{\rm{.012i}}} & {\rm{0}} & {\rm{0}}\\[5pt]{{\rm{0}}{\rm{.0069}} + {\rm{0}}{\rm{.0005i}}} & {\rm{0}} &{{\rm{0}}{\rm{.12}} - {\rm{0}}{\rm{.0011i}}} & {{\rm{ - 0}}{\rm{.084}} - {\rm{0}}{\rm{.0045i}}}\\[5pt]{{\rm{0}}{\rm{.052 - 0}}{\rm{.0028i}}} & {{\rm{1}}{\rm{.3}} - {\rm{0}}{\rm{.012i}}} & {{\rm{ - 0}}{\rm{.70}} + {\rm{0}}{\rm{.038i}}} & {{\rm{ - 0}}{\rm{.31}} + {\rm{0}}{\rm{.017i}}}\\[5pt]{{\rm{0}}{\rm{.16}} - {\rm{0}}{\rm{.0040i}}} & {\rm{0}} & {{\rm{ - 0}}{\rm{.41}} + {\rm{0}}{\rm{.016i}}} &{{\rm{ - 1}}{\rm{.3}} + {\rm{0}}{\rm{.031i}}}\end{array}} \right]\end{equation}

After obtained the FRF matrices $H\left( {18} \right)$ , $H\left( {36} \right)$ and $H\left( {54} \right)$ , the simulation of identifying the hub dynamic loads can be carried out by using the following four hub dynamic loads:

(23) \begin{align}{{F_x} = 3.0\sin \left( {2\pi \times 18t + \frac{\pi }{3}} \right) - 1.44\sin \left( {2\pi \times 36t + \frac{\pi }{6}} \right) - 0.46\sin \left( {2\pi \times 54t - \frac{\pi }{3}} \right)}\nonumber\\[8pt]{{F_z} = 10.0\sin \left( {2\pi \times 18t + \frac{\pi }{2}} \right) + 4.2\sin \left( {2\pi \times 36t + \frac{\pi }{3}} \right) + 1.83\sin \left( {2\pi \times 54t - \frac{\pi }{3}} \right)}\nonumber\\[8pt]{{M_y} = - 1.2\sin \left( {2\pi \times 18t + \pi } \right) + 0.41\sin \left( {2\pi \times 36t - \frac{\pi }{2}} \right) - 0.18\sin \left( {2\pi \times 54t - \frac{\pi }{6}} \right)}\nonumber\\[8pt]{{M_z} = 0.3\sin \left( {2\pi \times 18t - \frac{\pi }{6}} \right) + 0.09\sin \left( {2\pi \times 36t + \frac{\pi }{6}} \right) - 0.022\sin \left( {2\pi \times 54t + \frac{\pi }{2}} \right)}\end{align}

Table 4. Identified amplitudes and phases of hub dynamic loads at 18Hz

Table 5. Identified amplitudes and phases of hub dynamic loads at 36Hz

Under simultaneous excitation of the four dynamic loads by Equation (23), the vibration acceleration responses at the six measuring points on the frame structure can be obtained by using the Newmark-Beta method. Using the inverse of the FRF matrices by Equations (20), (21) and (22) to identify the dynamic loads by Equation (8), the identified results and errors of the dynamic loads at the three exciting frequencies can be obtained, respectively. The identified amplitudes and phases of the applied four dynamic loads and their errors are listed in Tables 4, 5 and 6. It can be seen from these tables that for the dynamic loads at the first three blade passage frequencies, the errors of the identified amplitudes are all within 2.1%, and the errors of identified phases are all within 5.2%, indicating that the identified dynamic loads have high accuracies, and the selected measuring points are reasonable.

Table 6. Identified amplitudes and phases of hub dynamic loads at 54Hz

3.3 Layout of measuring points by condition number of FRF matrix

In this section, the condition number of the FRF matrix is used to select the layout of the six measuring points (including position and direction). In order to indicate the influence of the layout of the six measuring points on the identification accuracy of the dynamic loads, the identification simulations of the four dynamic loads under different layout of the six measuring points were carried out. The amplitudes of the dynamic loads were set to 100N or 10Nm. 10,000 layouts including positions and directions of six measuring points were randomly generated. In each layout of six measuring points, the amplitudes of the four dynamic loads were identified separately. The average identification errors of the amplitudes of the four dynamic loads at the 10,000 layouts of six measuring points are shown in Fig. 5 in which the average identification errors are the average values of the identification errors at six measuring points in each layout. It can be seen from Fig. 5, the average identification errors at some layouts are very large, even reaching 5,000%. If the average error less than 5% is as qualified standard, then the qualification rate of the average identification errors in the randomly generated 10,000 layouts is only 37.48% according to the analysis of the average identification errors shown in Fig. 5. Therefore, it is necessary to reasonably select the layouts of six measuring points according to the condition number of the FRF matrix.

Figure 5. Amplitude identification errors of four dynamic loads at randomly generated 10,000 layouts of six measuring points.

To carry out the identification simulation of the four dynamic loads with reasonable layouts of six measuring points, 1,000 layouts of six measuring points were randomly generated by setting a range of the condition number of the FRF matrix. The qualification rate of the average identification errors and the maximum average identification error of the amplitude identification of the four dynamic loads at the 1,000 layouts of six measuring points in three ranges of the condition number of the FRF matrix are shown in Table 7.

Table 7. Amplitude identification results of four dynamic loads at randomly generated 1,000 layouts of six measuring points in three ranges of condition numbers

To clearly indicate the layout effect of six measuring points on identifying accuracy according to the condition number of the FRF matrix, four layouts of six measuring points were randomly generated in four ranges of the condition number of FRF matrix to identify the amplitudes of the four dynamic loads, respectively. The identification results and errors are listed in Table 8. It can be seen from Table 8 that with the decrease of the ranges of the condition number of FRF matrix, the amplitude identification error of the four dynamic loads decreases rapidly, indicating that the layout effect of six measuring points on identifying accuracy according to the condition number of the FRF matrix is very good. For the 4th range of the condition number, i.e. $290 \lt cond(H) \lt 300$ , the amplitude identification errors of the four dynamic loads are all less than 1%. Two layouts of six measuring points for the 1st and 2nd ranges of the condition number are shown in Figs. 6 and 7 respectively. It can be seen from Figs. 6 and 7 that the two layouts are quite different.

Table 8. Amplitude identification results and errors of four dynamic loads in four ranges of condition numbers

Figure 6. Layout of measuring points for condition number 2,900 < cond (H) < 3,000.

Figure 7. Layout of six measuring points for condition number 1,900 < cond (H) < 2,000.

4.1 Experimental system

The frame structure used for the experiment is shown in Fig. 1, and the schematic diagram of the experimental system for the identification of the hub dynamic loads is shown in Fig. 8. The experimental system included the PCB 208C05 force sensor, PCB 352C65 acceleration sensor, TMS320F28335 digital signal processor, HEV-50 shaker and so on. The vibration signals collected by the accelerometers were amplified to improve the signal-to-noise ratios of measured signals using the TMS320F28335 digital signal processor for analog-to-digital conversion and obtain high-quality vibration response acceleration data, and the amplification factor is set to 100. A dynamic force generated by the shaker at the excitation point of the frame structure was used to simulate the hub dynamic load, which produced the vibration of the frame structure. If the excitation point is on the neutral line of the frame structure, the dynamic load is a force in z direction. If the excitation point is not on the neutral line of the frame structure, the dynamic load will form a force in z direction and a torque in y direction. One excitation point and three measuring points in z direction on the frame structure shown in Fig. 7 were selected to carry out the identification experiment of the hub dynamic loads.

Figure 8. Schematic diagram of experimental system for identification of hub dynamic loads.

4.2 Experimental results

The excitation frequency in the experiment was set to 30Hz, which is different from the frequencies 18, 36 and 54Hz used in the simulation in order to further verify the feasibility of the proposed identification method in this paper for different frequencies. The FRF matrix between the excitation point and three measuring points was obtained under unit force excitation. The acceleration responses at three measuring points when excitation point was located at the neutral position and at 7cm from the neutral position were obtained, respectively. By subtracting the acceleration values at the same measuring point and at two excitation points, the acceleration response at the measuring point under the excitation of unit torque in y direction can be obtained. By analysing the response data at each measuring point, the FRF matrix between the excitation point and the three measuring points at 30Hz was obtained as follows:

(24) \begin{equation}H\left( {30} \right) = \left[ {\begin{array}{c@{\quad}c}{ - 2.98 - 2.14i} & {0.37 + 0.29i}\\[5pt]{ - 6.95 - 6.10i} & {4.44 + 0.054i}\\[5pt]{ - 4.66 - 4.04i} & {0.050 + 0.0099i}\end{array}} \right]\end{equation}

By analysing the FRF matrix (Equation (24)), the condition number of the matrix is $cond\left( {H\left( {30} \right)} \right) = 4.99$ , indicating that the selected position and direction of the three measuring points is a reasonable layout and the FRF matrix is well-conditioned. After obtaining the FRF matrix, the identification experiment of the dynamic loads can be carried out. The identification experiment of the dynamic loads included nine excitation conditions including a combination of three positions of excitation point deviating from the neutral line of the frame structure with three amplitudes of exciting force under each position of excitation point. The three positions of excitation point deviated 3.5, 5.25 and 7.0cm from the neutral line, respectively, and the three amplitudes of exciting force were 6.67, 8.33 and 10.00N, respectively. The experimental results and errors of the dynamic load identification in the nine excitation conditions are listed in Tables 9, 10 and 11, respectively. It can be seen from these tables that the amplitude identification errors of the dynamic force and torque in the experiment are less than 3%, indicating again that the identification of the hub dynamic loads containing force and torque based on the condition number of the FRF matrix measured on the fuselage structure can achieve an excellent effect.

It should be noted that although the frame structure used in this paper is simple comparing with a helicopter, the proposed method for identifying the hub dynamic loads in this paper is still suitable for a helicopter. The rotor hub dynamic loads are transmitted to the fuselage and cause vibration of the fuselage, and conversely can be identified by the vibration of the fuselage. Hence, the proposed method for identifying the hub dynamic loads in this paper is suitable for a helicopter by using Equation (8) to solve the hub dynamic loads as long as the measured FRF matrix of the measuring points on the fuselage meets the requirement of the condition number.

Table 9. Experimental results and errors of dynamic load identification when excitation point deviated 3.5cm from neutral line

Table 10. Experimental results and errors of dynamic load identification when excitation point deviated 5.25cm from neutral line

Table 11. Experimental results and errors of dynamic load identification when excitation point deviated 7.0 cm from neutral line