An accurate method for identifying hub dynamic loads by condition number of measured FRF matrix on helicopter fuselage

Abstract

It is a simple method to identify the hub dynamic loads of rotor by measuring the vibration responses on helicopter fuselage. However, the identification accuracy of the hub dynamic loads is related to the layout or placement of measuring points on the fuselage. The identification will be inaccurate if the layout of measuring points on the fuselage is unreasonable to result in the “ill-conditioned” frequency response function (FRF) matrix measured on the fuselage. In order to avoid the inaccurate identification due to the “ill-conditioned” measured FRF matrix, an accurate method for identifying the hub dynamic loads of rotor by vibration measurement on helicopter fuselage is proposed in this paper. In the proposed method, the reasonable layout of the measuring points on the fuselage for the “well-conditioned” measured FRF matrix can be obtained according to the condition number of the measured FRF matrix on the fuselage, and then the hub dynamic loads of rotor can be accurately identified. The simulation and experiment of the identification of the hub dynamic loads on a dynamically similar frame structure of a helicopter cockpit floor have verified the effectiveness and accuracy of the proposed method.

Nomenclature

$cond(H)$

condition number of the matrix $\left[ {H\left( \omega \right)} \right]$

[C]

damping coefficient matrix

${F_j}\left( \omega \right)$

exciting load at point j of the structure

${\left[ {F\left( \omega \right)} \right]_{{n_f} \times 1}}$

column vector of exciting loads

${F_x}$

lateral force

${F_y}$

longitudinal force

${F_z}$

vertical force

${H_{ij}}\left( \omega \right)$

FRF between points i and j

[H]

FRF matrix

${\left[ {H\left( \omega \right)} \right]^{ - 1}}$

inverse of the FRF matrix

[K]

stiffness matrix

${M_x}$

lateral moment

${M_y}$

longitudinal moment

${M_z}$

vertical moment

[M]

mass matrix

${n_f}$

number of exciting loads to be identified

${n_u}$

number of vibration measuring points on the structure

${N_b}$

number of blades

$\Delta t$

time increment

${U_i}\left( \omega \right)$

response generated at point i of the structure

${\left[ {U\left( \omega \right)} \right]_{{n_u} \times 1}}$

column vector of measured responses

x

lateral direction of the fuselage

y

longitudinal direction of the fuselage

z

vertical direction of the fuselage

Greek symbol

${\sigma _{\max }}(H)$

maximum singular value of the FRF matrix $\left[ {H\left( \omega \right)} \right]$

${\sigma _{\min }}(H)$

minimum singular value of the FRF matrix $\left[ {H\left( \omega \right)} \right]$

ω

frequency

Ω

rotor speed

Subscripts/Superscripts

H

complex conjugate transpose of matrix

+

generalised inverse of matrix

Abbreviations

FBG

fibre Bragg grating

FRF

frequency response function

1.0 Introduction

The dynamic loads of rotor hub are transmitted to helicopter fuselage through rotor shaft, causing strong vibration of the fuselage. Therefore, the accurate identification of the hub dynamic loads is of great significance for analysis and control of the fuselage vibration. The aerodynamic, inertia, elastic, damping and centrifugal forces of each rotor blade are synthesised to blade root, forming the root forces of each blade, then superimposing to the centre of the rotor hub to form the hub dynamic loads which are the three forces and three moments in three directions. At present, there are three methods for identification of the hub dynamic loads as follows.

The first method for identification of the hub dynamic loads is based on analysis of the rotor blade, i.e., analysing firstly the dynamic loads on each blade, then obtaining the dynamic loads at the blade root, and then obtaining the hub dynamic loads. This method requires establishment of complex aeroelastic dynamic model of the blade and its solving method. Because the aerodynamic environment and movement of the rotor blade are very complex during the flight of a helicopter, the established aeroelastic dynamic model of the rotor blade involves the elastic dynamics, unsteady aerodynamics, wake modeling and so on. At present, there are two kinds of the method for calculating the hub dynamic loads according to the different methods for calculating the aerodynamic loads of the blade and the wake. One is by using the comprehensive analysis models [1, 2] which are the aeroelastic dynamic models of the blade established by combining the multi-body dynamics and structural dynamics with the approximate aerodynamic loads and wake calculation methods. The typical examples are CAMRAD II [3], RCAS [4], DYMORE4.0 [5], UMARC2 [6], etc. The other is by adopting the CFD/CSD coupled aeroelastic dynamics. The typical representatives of commonly used CFD algorithms include OVERFLOW [7], TURNS [8] and FUN3D [9], etc. The prediction accuracy and capability of CFD/CSD coupled aeroelastic dynamics for the hub dynamic loads are significantly improved. However, under complex manoeuvering conditions and specific flight conditions, the CFD/CSD prediction accuracy for the harmonic components of the hub dynamic loads is still difficult to meet the requirements of fuselage vibration analysis [10].

The second method for identification of the hub dynamic loads is based on the vibration measurement on the rotor, i.e., measuring the dynamic strains at the blade root or the distributed dynamic strains on the blade or the dynamic strains on the rotor shaft, and then analysing the hub dynamic loads. The blade dynamic strains can be measured by using strain gauges or fibre Bragg grating (FBG) sensors. The measurement using strain gauges is the most widely used method for dynamic strain measurement on the blade and has been applied in lots of measurements of typical helicopters such as UH-60A [11]. In recent years, the FBG sensors have been used for dynamic strain measurements of the rotor blade because of the advantages such as anti-electromagnetic interference, high measurement accuracy, wide measurement range and quasi-distributed measurement [12, 13]. At present, there are some static calibration measurement methods and various load identification methods [14, 15] which can be used to identify the hub dynamic loads based on the strain data measured on the blade, but due to the complexity of practical problems, these methods adopt some specific assumptions, resulting in accuracy loss of the identified hub dynamic loads.

The third method for identification of the hub dynamic loads is based on the vibration measurement on the fuselage, i.e., measuring the vibration responses on the fuselage, and then analyzing the hub dynamic loads. In this method, the vibration acceleration responses on fuselage are firstly measured to obtain the FRF matrix between the exciting points and measuring points, then inverse the FRF matrix, and then obtain the hub dynamic loads. Giansante et al. [16] placed acceleration sensors on the fuselage and used the FRF matrix measured on the ground to identify the dynamic loads of AH-1G helicopter in flight. Fabunmi [17] studied the factors affecting the measurement accuracy of the helicopter dynamic loads, and experimentally showed that the modal parameters could be used to evaluate the reliability of dynamic loads measurement. Dobson et al. [18] studied the inverse problem of the FRF matrix used for identification of the dynamic loads, and proposed that the identification error could be reduced by increasing the number of measuring points. Callahan et al. [19] studied the identification of applied positions of dynamic loads based on the inverse of the FRF matrix. The results showed that a better identification effect could be achieved when the identification range included all actual applied positions of dynamic loads.

Among the above three methods, the third method is the easiest way to implement. However, in the analysis process, the FRF matrix obtained by the vibration measurements on the fuselage needs to be inversed. Hence the measured FRF matrix must be ‘well-conditioned’ for an accurate identification of the hub dynamic loads. However, the identification of the hub dynamic loads will be inaccurate if the measured FRF matrix is ‘ill-condition’. Whether the measured FRF matrix is well-conditioned or ill-conditioned depends on a reasonable or unreasonable layout of the measurement points on the fuselage since the position of the hub dynamic loads is determined. The reasonable layout of measurement points on the fuselage means that measurement points have a large vibration response, and the signal-to-noise ratios measured at measurement points are high, and then the measured FRF matrix will be well-conditioned. The unreasonable layout of measurement points on the fuselage means that the vibration responses at some measurement points are very small, and the signal-to-noise ratios measured at these measurement points are very low, and then the measured FRF matrix may be ill-conditioned. Therefore, it is necessary to obtain a reasonable layout of the measuring points on the fuselage to obtain a well-conditioned FRF matrix to improve the identification accuracy of the hub dynamic loads. At present, in the field of dynamic load identification, the problems of ill-conditioned FRF have been investigated [20, 21]. Jia et al. [20] proposed a weighted total least squares method to reduce the errors of the generally ill-posed random dynamic load identification. Alqam et al. [21] presented an approach for indirect identification of dynamic loads using the strain FRF, displacement FRF, along with the optimal locations.

It was difficult to obtain accurate identification of the hub dynamic loads by vibration measurement on the fuselage. The accuracy of identification of the hub dynamic loads was improved generally by adding measuring points on the fuselage [18]. According to matrix theory, whether a matrix is ill-conditioned depends on the condition number of the matrix. Hence, the condition number of the measured FRF matrix can be used as the basis of whether the layout of the measuring points on the fuselage is reasonable. Based on the matrix theory, an accurate method for identifying the hub dynamic loads by vibration measurement on helicopter fuselage is proposed in this paper. In the proposed method, the reasonable layout of the measuring points on the fuselage for the well-conditioned measured FRF matrix can be obtained according to the condition number of the measured FRF matrix on the fuselage, and then the hub dynamic loads of rotor can be accurately identified. The formula for calculating the condition number of the two norm FRF matrix is proposed in this paper. The simulations and experiments of the identification of the hub dynamic loads on a dynamically similar frame structure of Z-11 helicopter cockpit floor are carried out to verify the effectiveness and accuracy of the proposed method.

2.0 Accurate method for identifying hub dynamic loads based on condition number of measured FRF on fuselage

Helicopter fuselage vibration is mainly caused by the hub dynamic loads. The fuselage vibration frequencies, i.e., the frequencies of the hub dynamic loads are the blade passage frequencies $k{N_b}\Omega $ , where $k = 1,2,3, \ldots $ , ${N_b}$ is the number of blades, and Ω is rotor speed. Therefore, the fuselage vibration is the deterministic harmonic vibration.

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.0 Simulation of identifying hub dynamic loads on fuselage model structure

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.0 Experiment of identifying hub dynamic loads on fuselage model structure

In order to further verify the identification effectiveness and accuracy of the hub dynamic loads by using the condition number of the FRF matrix measured on fuselage, the experiment of identifying the hub dynamic loads is carried out in this section.

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

5.0 Conclusions

A simple and accurate method for identifying the hub dynamic loads based on the condition number of the FRF matrix obtained by the vibration measurements on helicopter fuselage has been proposed in this paper. The proposed method can effectively select the reasonable layouts of the measuring points on the fuselage, which can avoid the ill-conditioned measured FRF matrix, hence can significantly improve the identification accuracy of the hub dynamic loads. The specific conclusions of the identification simulation and experiment by using the proposed method in this paper are as follows:

(1) The identification simulations of the hub dynamic loads on a dynamically similar frame structure of Z-11 helicopter fuselage cockpit floor have small identification errors of the amplitudes and phases of four dynamic loads within 2.1% and 5.2%, respectively, under the first three blade passage frequencies of Z-11 helicopter. As the condition numbers of the FRF matrix decrease, the amplitude identification errors of the four dynamic loads decrease rapidly, indicating that the layout efficiency of the measuring points on the fuselage based on the condition number of the FRF matrix is very high.

(2) The identification experiments of the hub dynamic loads on a dynamically similar frame structure of Z-11 helicopter cockpit floor have small identification errors of the amplitudes of the dynamic force and torque in nine exciting conditions less than 3%, experimentally verifying the identification accuracy of the hub dynamic loads containing force and torque based on the condition number of the FRF matrix measured on the fuselage.