Model Parameter Estimation Method on Multichannel Structural Response from Turbulence Flutter Test

Abstract

Flutter design is important for the design of new aircraft types, for which flutter tests are an important verification measure. Atmospheric turbulence excitation is a common form of excitation in flutter flight tests. The modal parameter estimation of the turbulence response is a key aspect to ensure accurate data processing of flutter test results. Atmospheric turbulence excitation acts on the structural system, and the turbulence response thus simultaneously contains both the randomness of the excitation signal and the determinism of the structure system. In view of the turbulence response characteristics, this paper addresses the incoherence of atmospheric turbulence excitation and the orthogonality of the frequency domain from a multichannel response. The turbulence response is used to perform modal parameter identification in the frequency domain. The power spectral density matrix can be calculated from the multichannel turbulence response using the periodogram method. Singular value decomposition is then performed on the power spectral density matrix at each spectral pin based on the orthogonality of the frequency domain. The maximum singular value of each spectral pin forms a curve over the entire frequency band, which is the autopower spectral density function of the system, the system is directly identified at the frequency domain using the polynomial fitting in the frequency domain, and the modal parameters (frequency, damping ratio) are calculated according to the fitted transfer function. This paper verifies the theoretical feasibility of the proposed method using simulation data. The engineering applicability is verified based on the turbulence response from the flutter flight test of a certain aircraft type.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Flutter design is an important requirement in the overall design of new aircraft types. A sufficient flutter envelope prediction guarantees the safe flight of the aircraft. Flutter flight tests are thus a key link to probe the flutter design to verify unstable aeroelastic phenomena. The test signals collected in flight tests are an important basis for calculating the natural frequency and damping ratio of the dynamic aeroelastic structure at different velocities. However, factors such as experimental measurement noise, lack of excitation signal and excitation noise, and a limited number of data acquisition sensors can produce noise-contaminated data, which can lead to spurious mode problems in the modal parameter estimation.

Flutter flight tests and a variety of excitation methods have been developed alongside advances in aircraft test technology. According to the principle of system identification, the excitation for a structural system includes the pulse excitation, sweep frequency signal excitation, and atmospheric turbulence excitation in the flutter flight tests.

The impulse response caused by the pulse excitation is the most ideal signal for structural system identification. A small rocket excitation and rudder vibration are often used in flutter flight test to approximate the impulse excitation signals. The modal parameters can be determined by collecting the structure response. Thehabey [1] used the matrix pencil (MP) method to the dynamic instability analysis of aeroelastic systems. Kiviaho et al. [2] performed a study based on the complex exponential method and MP method to identify the modal parameters of the flutter test signal to predict the flutter boundary.

Sweep frequency excitations are a way to produce narrow-band excitation signal excitations near the natural frequency of the structure. This method is conducive to obtaining a high-quality structural response, and the excitation signal can be measured because the sweep frequency excitation is generated by the excitation source installed on the aircraft. Richardson and Formenti [3] used the rational fraction polynomial for frequency response function (FRF) identification from the perspective of system identification and applied the FRF to calculate the structural system modal parameters. In terms of the swept frequency excitation in flutter tests, the frequency polynomial fitting method can be directly used for system identification and modal parameter calculations because the excitation and response signals can be simultaneously measured [4].

The pulse excitation in flutter tests is often realized by a small rocket excitation or rudder operation. The risk of the flutter test itself is thus an important challenge for producing excitation action in terms of test safety. Sweep frequency excitations are a kind of narrow-band excitation signal, but the high-frequency response of the aircraft structure is often difficult to obtain due to the bandwidth of the excitation source. Atmospheric turbulence is a kind of random excitation. The frequency domain band of the excitation signal is sufficiently wide, while the energy of each frequency domain of the excitation signal is dispersed due to the bandwidth. This feature is mainly manifested in the turbulence response, in which the response signal contains both the information of the structural system and the randomness due to the random excitation.

Because turbulence excitation is not measurable, the system identification problem for the turbulence response is referred to as output-only modal parameter estimation. The modal parameter estimation of flutter response signals mainly includes the free decay response calculation method based on the random decrement technique (RDT), the stochastic subspace identification (SSI) method based on state space modeling, and the transfer function modeling method of the frequency domain decomposition (FDD). Wen et al. [5] used the RDT for modal identification of structures with close modes. Oymak and Ozay [6] revisited the Ho-Kalman-based system identification method and analyzed the robustness and performance of the finite sample. Scionti et al. [7] introduced the output-only response modal parameter identification method of SSI balanced realization and introduced the stochastic model based on the state space model.

Modal parameter estimation methods of the turbulence response signals based on transfer function modeling mainly include two types. (1) Frequency-time domain estimation methods are mainly based on FDD, which mostly considers the multichannel turbulence response to obtain the power spectral density (PSD) matrix. The auto power spectrum of a single degree-of-freedom (SDOF) system can be calculated using the FDD method to obtain the autocorrelation function of the SDOF system using the inverse FFT. The model parameters are then calculated in the time domain [8, 9]. (2) The least-square complex frequency-domain (LSCF) method is used to model the parameter estimation of the aircraft [10] and to study the online monitoring of aircraft modal parameter analysis using only known outputs [11].

According to the transfer function modeling analysis of the turbulence response, the FDD method must perform an inverse FFT for the auto PSD function of a SDOF system to identify the modal parameters. The LSCF method focuses on the turbulence of multichannels, and the turbulence response is based on the LSCF estimator to model the common-denominator transfer function. The modal parameters are thus calculated using the identified transfer function, and a stability diagram is used to estimate the stable modal parameters.

This paper proposes to calculate the PSD function of the multichannel turbulence response based on the FDD method. The singular value decomposition (SVD) for the PSD function matrix is performed at each spectral pin. The maximum singular value curve is mathematically given as the product of the system transfer function and its conjugate transpose. The maximum singular value curve can therefore be fitted in the frequency domain using the rational fraction polynomial and calculated the modal parameters. The contributions of the method can be summarized as follows.

(1) The influence of random noise on the spectrum estimation is reduced using multichannel power spectrum analysis and SVD based on the randomness of the turbulence response signal in the flutter test

(2) The system identification and modal parameter estimation are performed by frequency domain curve fitting based on the relationship between the maximum singular value curve and the system transfer function

(3) The FRF is fitted using the rational fraction polynomial method at the analysis frequency band to improve the performance of the modal parameter estimation

(4) The method is verified based on the turbulence response signal of the flutter flight test, which verifies the feasibility and engineering applicability of multimodal frequency domain fitting

2. Related Methods

2.1. Singular Value and Transfer Function

Under ideal conditions, the following relationship exists between the output $yt$, input $xt$, and the impulse response $ht$ of the system [12]: $$\begin{array}{}\text{(1)}& yt={\int }_0^{\infty }h\tau xt-\tau d\tau ,\end{array}$$where if $\tau <0$, $h\tau =0$ is a physically realizable system, and then $ytyt+\tau$ can be written as follows: $$\begin{array}{}\text{(2)}& ytyt+\tau ={\int }_0^{\infty }{\int }_0^{\infty }h\alpha h\beta xt-\beta xt+\tau -\alpha d\alpha d\beta .\end{array}$$

Eq. (2) can be written in the form of an input and output autocorrelation: $$\begin{array}{}\text{(3)}& R_{yy}\tau ={\int }_0^{\infty }{\int }_0^{\infty }h\alpha h\beta R_{xx}\tau +\beta -\alpha d\alpha d\beta .\end{array}$$

Similarly, $xtyt+\tau$ also can be written as follows: $$\begin{array}{}\text{(4)}& xtyt+\tau ={\int }_0^{\infty }h\alpha xtxt+\tau -\alpha d\alpha .\end{array}$$

Furthermore, Eq. (4) can be expressed as follows: $$\begin{array}{}\text{(5)}& R_{xy}\tau ={\int }_0^{\infty }h\alpha R_{xx}\tau -\alpha d\alpha .\end{array}$$

Eqs. (3) and (5) are a convolution form, as in Eq. (1), and the double-sideband spectral densify functions $S_{xx}j\omega$, $S_{yy}j\omega$, and $S_{xy}j\omega$ can be calculated by Fourier transform: $$\begin{array}{}\text{(6)}& S_{yy}j\omega ={Hj\omega }^2{S}_{xx}j\omega ,\text{(7)}& S_{xy}j\omega =Hj\omega S_{xx}j\omega .\end{array}$$

Because Eqs. (6) and (7) are double-sideband spectral density functions, the frequency ranges from positive to negative. Eq. (6) is a real-valued relation that contains only the gain factor $Hj\omega$ of the system. In Eq. (7), $Hj\omega$ contains both the amplitude spectrum $Hj\omega$ and phase spectrum $\varphi j\omega$. Equation (7) thus becomes the correlation between the input and output cross-spectrum. The input-output relationship of a system in an ideal state has been proved, and the system is a linear time-invariant system. The relationship between the input/output autocorrelation function and transfer function can be more intuitively understood by means of the frequency domain conversion. The relationship between the double-sideband autocorrelation and transfer function can therefore be determined using the one-sideband spectral density function as follows: $$\begin{array}{}\text{(8)}& G_{yy}j\omega ={Hj\omega }^2{G}_{xx}j\omega ,\text{(9)}& G_{xy}j\omega =Hj\omega G_{xx}j\omega .\end{array}$$

Compared with Eq. (1), $$\begin{array}{}\text{(10)}& Yj\omega =Hj\omega Xj\omega ,\end{array}$$where $Xj\omega$ and $Yj\omega$ are the Fourier transform of $xt$ and $yt$, respectively. The complex frequency domain of Eq. (10) exists as follows: $$\begin{array}{c}\text{(11)}& Y^{\ast }j\omega =H^{\ast }j\omega X^{\ast }j\omega ,{Yj\omega }^2={Hj\omega }^2{Xj\omega }^2,\\ X^{\ast }j\omega Yj\omega =Hj\omega {Xj\omega }^2.\end{array}$$

The complex conjugate of Eq. (9) is thus as follows: $$\begin{array}{}\text{(12)}& G_{xy}^{\ast }j\omega =G_{yx}j\omega =H^{\ast }j\omega G_{xx}j\omega ,\end{array}$$where $$\begin{array}{}\text{(13)}& G_{yx}j\omega =G_{yx}j\omega e^{j{\theta }_{xy}j\omega },\text{(14)}& H^{\ast }j\omega =Hj\omega e^{j\varphi j\omega }.\end{array}$$

For the system FRF in the complex domain, Eqs. (8) and (13) can be further written as $$\begin{array}{}\text{(15)}& G_{yy}j\omega =Hj\omega H^{\ast }j\omega G_{xx}j\omega .\end{array}$$

According to the relationship between the input and output spectral density functions and the FRF of the system, the following can be written: $$\begin{array}{}\text{(16)}& G_{yy}j\omega =Hj\omega G_{xx}j\omega H{j\omega }^H,\end{array}$$where $G_{xx}j\omega$ is a $r\times r$ PSD function matrix of the input signal, $r$ is the channel number of the input signal, $G_{yy}j\omega$ is an $m\times m$ PSD function matrix of the output signal, $m$ is the channel number of the output signal, and $Hj\omega$ is an $m\times r$ FRF matrix. After decomposing $Hj\omega$ using the residue method, it can be expressed as follows: $$\begin{array}{}\text{(17)}& Hj\omega =\sum _{k=1}^n\frac{R_k}{j\omega -{\lambda }_k}+\frac{{\overline{R}}_k}{j\omega -{\overline{\lambda }}_k},\end{array}$$where $n$ is the model number, ${\lambda }_k$ is the pole, and $R_k$ is the residual: $$\begin{array}{}\text{(18)}& R_k={\varphi }_k{\gamma }_k^T,\end{array}$$where ${\varphi }_k$ is the model shape vector, and ${\gamma }_k$ is the model participation vector according the irrelevance of the multichannel input excitation signal. The PSD matrix of the excitation is a constant matrix, $G_{xx}j\omega =C$. Equation (16) can therefore be written as follows: $$\begin{array}{}\text{(19)}& G_{yy}j\omega =\sum _{k=1}^n\sum _{s=1}^n\frac{R_k}{j\omega -{\lambda }_k}+\frac{{\overline{R}}_k}{j\omega -{\overline{\lambda }}_k}C{\frac{R_s}{j\omega -{\lambda }_s}+\frac{{\overline{R}}_s}{j\omega -{\overline{\lambda }}_s}}^H,\end{array}$$where Eq. (19) can be obtained by the Heaviside partial fraction expansion theory: $$\begin{array}{}\text{(20)}& G_{yy}j\omega =\sum _{k=1}^n\frac{A_k}{j\omega -{\lambda }_k}+\frac{\overline{A_k}}{j\omega -{\overline{\lambda }}_k}+\frac{B_k}{-j\omega -{\lambda }_k}+\frac{\overline{B_k}}{-j\omega -{\overline{\lambda }}_k},\end{array}$$where $A_k$ represents the residual matrix of the $k$-^th mode of the response signal correlation PSD matrix, and $A_k$ is a Hermitian matrix with dimensions of $m\times m$. $$\begin{array}{}\text{(21)}& A_k=R_{k}C\sum _{k=1}^n\frac{{\overline{R_s}}^T}{-{\lambda }_k-{\overline{\lambda }}_s}+\frac{R_s^T}{-{\lambda }_k-{\lambda }_s}.\end{array}$$

The residual can be written as follows: $$\begin{array}{}\text{(22)}& A_k=\frac{R_{k}C{\overline{R_k}}^T}{2{\sigma }_k},\end{array}$$where ${\sigma }_k$ is the opposite of the real part of the pole, where ${\lambda }_k=-{\sigma }_k+j{\mu }_k$. When the damping is a small value in the $k$-^th mode, the residual part of the mode vector can be expressed as follows: $$\begin{array}{}\text{(23)}& A_k\propto R_{k}C\overline{R_k}=\varphi {{}_k\gamma }_k^{T}C\gamma {{}_k\varphi }_k^T=d_k\varphi {{}_k\varphi }_k^T,\end{array}$$where $d_k$ is a constant. In the frequency domain, only a part of the key mode frequency domain has an important influence on the FRF of the system, such as a typical first- or second-order mode. The frequency is marked as $\text{sub}\omega$. Therefore, in the case of small damping, the response spectral density function of the system can be expressed as follows: $$\begin{array}{}\text{(24)}& G_{yy}j\omega =\sum _{k\in \text{sub}\omega }\frac{d_k{\varphi }_k{\varphi }_k^T}{j\omega -{\lambda }_k}+\frac{{\overline{d}}_k{\overline{\varphi }}_k{\overline{\varphi }}_k^T}{j\omega -{\overline{\lambda }}_k}\end{array}$$where $d_k$ is a constant, ${\varphi }_k$ is the mode shape, and ${\lambda }_k$ represents the pole. For the response signal of multiple measurement points, the PSD function ${\widehat{G}}_{yy}j\omega$ of the response signal is subjected to SVD at a discrete spectral pin in the following form: $$\begin{array}{}\text{(25)}& {\widehat{G}}_{yy}j\omega =U_i{S}_i{U}_i^T,\end{array}$$where $U_i=u_{i1},u_{i2},\cdots ,u_{im}$ is a unitary matrix of singular value vector $u_{ij}$, and $S_i$ is a diagonal matrix of singular values $s_{ij}$ for the $k$-^th mode. The corresponding frequency domain amplitude is reflected in the singular value. Compared with Eqs. (24) and (25), the first singular value vector $u_{i1}$ represents the mode shape, namely, $$\begin{array}{}\text{(26)}& \widehat{\varphi }=u_{i1}.\end{array}$$

This indicates that for a randomly excited system, the power spectrum of the multichannel excited signal is a flat spectrum, and the correlation power spectrum density function matrix is a diagonal matrix. A comparison of Eqs. (24) and (25) shows that the maximum singular value of the PSD function matrix ${\widehat{G}}_{yy}j\omega$ of the response at each spectral pin can be expressed as the auto PSD function of the corresponding system.

2.2. Rational Fraction Polynomial

The rational fraction polynomial is an orthogonal polynomial method used to perform a frequency polynomial fitting of the maximum singular value curve. The model parameters can be calculated from the estimated FRF based on the rational fractional polynomial.

The FRF for a linear time-invariant system can be expressed as follows: $$\begin{array}{}\text{(27)}& Hs=\frac{a_m{s}^m+\cdots +a_2{s}^2+a_{1}s+a_0}{b_n{s}^n+\cdots +b_2{s}^2+b_{1}s+b_0},\end{array}$$where $m<n$. Dividing the numerator and denominator of Eq. (27) by $b_n$ gives $$\begin{array}{}\text{(28)}& \begin{array}{c}{p}_{m}s=s^m,\cdots ,p_{1}s=s,p_{0}s=1,\\ q_{n}s=s^n,\cdots ,q_{1}s=s,q_{0}s=1.\end{array}\end{array}$$

This leads to the following: $$\begin{array}{}\text{(29)}& Hs=\frac{{\sum }_{k=0}^m{c}_k{p}_{k}s}{{\sum }_{k=0}^n{d}_k{q}_{k}s},\text{where}{d}_n=1.\end{array}$$

In practical engineering, it is often impossible to measure the FRF over the entire $S$-plane; thus, only the data that cover the entire frequency axis are considered. These data are generally known as the FRF, i.e., $s=j\omega$. In this case, Eq. (29) can be further expressed as follows: $$\begin{array}{}\text{(30)}& Hj\omega =\frac{{\sum }_{k=0}^m{c}_k{p}_{k}j\omega }{{\sum }_{k=0}^n{d}_k{q}_{k}j\omega },\end{array}$$where the coefficients $c_k$ and $d_k$ are the quantities that need to be calculated. To use the orthogonal polynomial to fit the FRF, it is necessary to simplify the calculation using the characteristics of the orthogonal polynomial and the FRF. In general, only the FRF calculation in the positive frequency range is considered; however, the FRF in the negative frequency range also exists, which implies that the FRF is a conjugate symmetric function about the origin of the frequency axis. The concept of negative frequency is therefore introduced in the calculation process, in which $\omega ={\omega }_{-L}\cdots {\omega }_{-1}{\omega }_1\cdots {\omega }_L$ over a total of 2 L points. If we let ${\omega }_{-i}=-{\omega }_i$, then $Hj{\omega }_{-i}=H‐j{\omega }_i=H^{\ast }j{\omega }_i$, where $i=1,2,\cdots ,L$, and the measured value of the FRF $\widehat{H}j{\omega }_i$ exists according to the following: $$\begin{array}{}\text{(31)}& \widehat{H}j{\omega }_{-i}=\widehat{H}-j{\omega }_i=\widehat{H^{\ast }}j{\omega }_{i}i=1,2,\cdots ,L.\end{array}$$

The error between the theoretical model of the FRF (i.e., a rational fraction of Eq. (30)) and the actual measurement result can be expressed as follows: $$\begin{array}{}\text{(32)}& {\ell }_i=Hj{\omega }_i-\widehat{H}j{\omega }_i=\frac{{\sum }_{k=0}^m{c}_k{p}_{k}j{\omega }_i}{{\sum }_{k=0}^n{d}_k{q}_{k}j{\omega }_i}-\widehat{H}j{\omega }_i.\end{array}$$

To linearize the coefficients $c_k$ and $d_k$, Eq. (32) can be rewritten as follows: $$\begin{array}{}\text{(33)}& e_i={\ell }_i\sum _{k=0}^n{d}_k{q}_{k}j{\omega }_i=\sum _{k=0}^m{c}_k{p}_{k}j\omega -\widehat{H}j{\omega }_i\sum _{k=0}^{n-1}{d}_k{q}_{k}j{\omega }_i+q_{n}j{\omega }_i.\end{array}$$

The error is defined by Eq. (34): $$\begin{array}{}\text{(34)}& J=\sum _{i=-L}^L{e}_i^{\ast }{e}_i=E^{H}E,\end{array}$$where $E={e_{-L},\cdots ,e_{-1},e_1,\cdots ,e_L}^T$, and $E^H$ is a conjugate transpose of $E$. Considering Eqs. (33) and (34), the vector of the error is defined as follows: $$\begin{array}{}\text{(35)}& E=PC-QD-W,\end{array}$$where $$\begin{array}{c}\text{(36)}& P=\begin{array}{c}{P}^{-}\\ P^{+}\end{array}=\begin{array}{cccc}{p}_{0}j{\omega }_{-L}& p_{1}j{\omega }_{-L}& \cdots & p_{m}j{\omega }_{-L}\\ ⋮& ⋮& \cdots & ⋮\\ p_{0}j{\omega }_{-1}& p_{1}j{\omega }_{-1}& \cdots & p_{m}j{\omega }_{-1}\\ p_{0}j{\omega }_1& p_{1}j{\omega }_1& \cdots & p_{m}j{\omega }_1\\ ⋮& ⋮& \ddots & ⋮\\ p_{0}j{\omega }_L& p_{1}j{\omega }_L& \cdots & p_{m}j{\omega }_L\end{array},Q=\begin{array}{c}{Q}^{-}\\ Q^{+}\end{array}=\begin{array}{cccc}\widehat{H}j{\omega }_{-L}{q}_{0}j{\omega }_{-L}& \widehat{H}j{\omega }_{-L}{q}_{1}j{\omega }_{-L}& \cdots & \widehat{H}j{\omega }_{-L}{q}_{n-1}j{\omega }_{-L}\\ ⋮& ⋮& \cdots & ⋮\\ \widehat{H}j{\omega }_{-1}{q}_{0}j{\omega }_{-1}& \widehat{H}j{\omega }_{-1}{q}_{1}j{\omega }_{-1}& \cdots & \widehat{H}j{\omega }_{-1}{q}_{n-1}j{\omega }_{-1}\\ \widehat{H}j{\omega }_1{q}_{0}j{\omega }_1& \widehat{H}j{\omega }_1{q}_{1}j{\omega }_1& \cdots & \widehat{H}j{\omega }_1{q}_{n-1}j{\omega }_1\\ ⋮& ⋮& \ddots & ⋮\\ \widehat{H}j{\omega }_L{q}_{0}j{\omega }_L& \widehat{H}j{\omega }_L{q}_{1}j{\omega }_L& \cdots & \widehat{H}j{\omega }_L{q}_{n-1}j{\omega }_L\end{array},\\ C={c_0,c_1,\cdots ,c_m}^T,\\ D={d_0,d_1,\cdots ,d_{n-1}}^T,\\ W=\begin{array}{c}{W}^{-}\\ W^{+}\end{array}={\stackrel{\wedge }{H}j{\omega }_{-L}{q}_{n}j{\omega }_{-L},\cdots ,\stackrel{\wedge }{H}j{\omega }_{-1}{q}_{n}j{\omega }_{-1},\stackrel{\wedge }{H}j{\omega }_1{q}_{n}j{\omega }_1,\stackrel{\wedge }{H}j{\omega }_L{q}_{n}j{\omega }_L}^T.\end{array}$$

To determine the minimum error value $J$, the partial derivative of $J$ with respect to the vectors $C$ and $D$ can be set to zero, as shown here: $$\begin{array}{}\text{(37)}& \frac{\partial J}{\partial C}=0,\frac{\partial J}{\partial D}=0.\end{array}$$

According to Eqs. (35) and (37), the matrix relationship in Eq. (38) can be defined as follows: $$\begin{array}{}\text{(38)}& \begin{array}{cc}Y& X\\ X^T& Z\end{array}\begin{array}{c}C\\ D\end{array}=\begin{array}{c}G\\ F\end{array},\end{array}$$where $$\begin{array}{}\text{(39)}& \begin{array}{c}Y=\frac{1}{2}{P}^{H}P+P^T{P}^{\ast },\\ Z=\frac{1}{2}{Q}^{H}Q+Q^T{Q}^{\ast },\\ G=\mathrm{Re}{P}^{H}W,\\ F=\mathrm{Re}{Q}^{H}W.\end{array}\end{array}$$

In Eq. (38), if $pj\omega$ and $qj\omega$ meet the following conditions, $$\begin{array}{c}\text{(40)}& \sum _{i=-L}^L{p}_s^{\ast }j{\omega }_i{p}_{t}j{\omega }_i=\sum _{i=-L}^L{p}_{s}j{\omega }_i{p}_t^{\ast }j{\omega }_i=\begin{array}{cc}1\hfill & s=t\hfill \\ 0\hfill & s\ne t\hfill \end{array}s,t=0,1,\cdots ,m,\sum _{i=-L}^L{\stackrel{\wedge }{H}j{\omega }_i}^2{q}_s^{\ast }j{\omega }_i{q}_{t}j{\omega }_i=\sum _{i=-L}^L{\stackrel{\wedge }{H}j{\omega }_i}^2{q}_{s}j{\omega }_i{q}_t^{\ast }j{\omega }_i=\begin{array}{cc}1\hfill & s=t\hfill \\ 0\hfill & s\ne t\hfill \end{array}s,t=0,1,\cdots ,n-1.\end{array}$$Then, the $Y$ and $Z$ matrices in Eq. (38) create the unit matrix, the final form of which is given as follows: $$\begin{array}{}\text{(41)}& \begin{array}{cc}{I}_1& X\\ X^T& I_2\end{array}\begin{array}{c}C\\ D\end{array}=\begin{array}{c}G\\ F\end{array}.\end{array}$$

The $C$ and $D$ matrices can be obtained by expanding Eq. (41) to the following: $$\begin{array}{}\text{(42)}& \begin{array}{c}I-X^{T}XD=F-X^{T}G,\\ C=G-XD.\end{array}\end{array}$$

Satisfying the condition of Eq. (40) makes $p_{k}j\omega$ and $q_{k}j\omega$ an orthogonal polynomial, which in turn makes the diagonal of the coefficient matrix in Eq. (38) an identity matrix. Equation (42) can then be used to separately solve the coefficients $C$ and $D$. Considering the derivation described, the Forsythe complex orthogonal polynomial can be used to reconstruct $p_{k}j\omega$ and $q_{k}j\omega$ [13].

3. Simulation Signal Verification

The simulation signal was used for the model parameter identification to verify the proposed method. The impulse response of the structural system is described as follows: $$\begin{array}{}\text{(43)}& ht=\sum _{i=1}^N{A}_i{e}^{-{\xi }_i{\omega }_{i}t}\sin {\omega }_{i}t+{\phi }_i,\end{array}$$where $N$ is the number of modes in the system, $A_i$ represents the amplitude of the $i$-^th model, ${\xi }_i$ represents the damping of the $i$-^th mode, ${\omega }_i$ represents the frequency of the $i$-^th mode, and ${\phi }_i$ represents the phase corresponding to the $i$-^th mode, which is used to simulate the measured signal from a different sensor on the same structural component. The turbulence response of the system is then determined according to the following: $$\begin{array}{}\text{(44)}& yt=ht\ast xt,\end{array}$$where $yt$ represents the turbulence response, $ht$ is the impulse response of the structural system, and $xt$ is the Gaussian white noise in the simulation.

Based on the simulation signal calculation method, the simulation signal parameter setting is carried out according to the principle of the second-order modal coupling, which is the cause of the flutter generation, and the modal parameter identification is carried out using the proposed method. The specific parameters of the three simulation signal settings are shown in Table 1.

Table 1

Parameters of the simulation signal.

Simulation signal	Amplitude	Frequency	Damping	Phase
Signal #1 (two-order)	0.8	5.65	0.0214	0.68π
	1	13.75	0.01505	0.33π
Signal #2 (two-order)	0.9	7	0.01514	0.51π
	1.5	11	0.00805	0.69π
Signal #3 (three-order)	0.8	3.75	0.03205	0.23π
	1.1	7.5	0.00305	0.51π
	0.5	11.25	0.02514	0.9π

The two-order modal frequencies in the simulation signal (signal #1) are 5.65 Hz and 13.75 Hz with respective amplitudes of 0.8 and 1, corresponding damping of 0.0214 and 0.01505, and phases of 0.68π and 0.33π. The turbulence response of the different measurement points is shown in Figures 1 and 2 that show the function of the corresponding turbulence response.

[figure(s) omitted; refer to PDF]

According to the theoretical analysis in Section 2, the PSD of the simulated signal of each channel can be calculated, and the maximum singular value curve of the entire frequency band is obtained (Figure 3). The model parameters can be estimated based on the rational fraction polynomials based on the maximum singular value curve. Figure 4 (top) shows the waveform of the PSD function of the turbulence response signal from a certain channel in the analysis frequency band, and Figure 4 (bottom) compares the maximum singular value curve and results of the system identification.

[figure(s) omitted; refer to PDF]

On the basis of the two-order modal simulation signals (signal #1) of 5 Hz and 13 Hz, the two-order models are further set to 7.5 Hz and 11.25 Hz to verify the two-order modal coupling in the flutter process.

Figure 5 shows the turbulence response signal of the two-order mode (signal #2), including the 7 Hz and 11 Hz frequency model, which, respectively, have amplitudes of 0.9 and 1.5, damping of 0.01514 and 0.00805, and phases of 0.51π and 0.69π. Figure 6 shows the PSD function of the time domain turbulence response. Figure 7 is the maximum singular value curve of the three-channel turbulence response. Figure 8 (top) shows the PSD of the turbulence response signal, and Figure 8 (bottom) compares the maximum singular value curve and the identification result of the rational fraction polynomial. A comparison of the top and bottom graphs in Figure 8 further illustrates the relationship between the maximum singular value curve and PSD of the turbulence response signal. The system identification of the turbulence response with the maximum singular value curve can thus be realized.

[figure(s) omitted; refer to PDF]

The two settings of the simulation signals verify the system identification method proposed in this paper from the perspective of the flutter generation principle. However, the test speed is often far from the true flutter speed considering test safety concerns. In actual engineering, there is a real need for the modal parameter identification of the turbulence response using coexisting multimodels and dispersed energy in each modal.

The method is therefore further verified using three models (signal #3). Figure 9 shows the turbulence response of the three channels including the three-mode frequencies of 3.75 Hz, 7.5 Hz, and 11.25 Hz, which, respectively, have amplitudes of 0.8, 1.1, and 0.5, damping of 0.03205, 0.00305, and 0.02514, and phases of 0.23π, 0.51π, and 0.9π. Figure 10 is the PSD corresponding to the turbulence response signal of the three channels, and Figure 11 is the maximum singular value curve of the PSD matrix calculated from the three-channel turbulence response. Figure 12 (top) shows one of the PSD of the three-mode turbulence responses, and Figure 12 (bottom) compares the maximum singular value curve of the three-mode turbulence response signal and the system identification result of the rational fraction polynomial. The final curve fitting results demonstrate that the proposed method can also identify the modal parameters for the turbulence response signals containing multimodes.

[figure(s) omitted; refer to PDF]

The system identification based on the rational fraction polynomial method is carried out on the basis of the signal waveform and system identification curve, and the modal parameter calculation is performed based on the identified system transfer function. Table 2 shows the parameter identification results of the three sets of simulation signals. The findings indicate that the error between the identification results and true values is relatively small.

Table 2

Model parameter identification results of the simulation signal.

	True value		Our method
	Frequency (Hz)	Damping	Frequency (Hz)	Damping
Signal #1	5.65	0.0214	5.6367	0.0185
	13.75	0.01505	13.736	0.0147
Signal #2	7.5	0.01514	7.4385	0.0135
	11.25	0.00805	11.2504	0.0089
Signal #3	3.75	0.03205	3.754	0.0371
	7.5	0.00305	7.4975	0.0036
	11.25	0.02514	11.2143	0.0260

4. Physical Test Verification

4.1. Model Parameter Estimation

The response signal excited by the atmosphere turbulence of the flutter flight test of a certain type of aircraft is used as the physical test signal to verify the proposed method. The acceleration response and corresponding PSD (Figures 13–16) locate the front and rear edges of the horizontal tail, as well as the front and rear edges of its wing tip.

[figure(s) omitted; refer to PDF]

The maximum singular value curve of the corresponding signal is calculated using the proposed method based on the turbulence response of the multiple measurement channels. The system identification is performed based on the rational fraction polynomial according to the maximum singular value curve, and the modal parameters are calculated. Figure 17 shows the maximum singular value curve of the four-channel response signals of the horizontal tail. The singular value curve clearly contains two modes at 11 Hz and 25 Hz. The system is identified using the rational fraction polynomial method based on the maximum singular value curve. The blue curve in Figure 18 is the maximum singular value curve, and the red curve is the system identification result. Due to the noise problem of the maximum singular value, it is not possible to consider only the two modes of interest in the order setting of the system identification. It is thus often necessary in engineering applications to increase the order to achieve better system identification and fit the modes of interest.

[figure(s) omitted; refer to PDF]

The turbulence response of the front and rear edges of the wing and front and rear edges of the wing tip is used to verify the proposed method. Figures 19–22 show a time series and the PSD function of the turbulence signal from the four different test channels of the wing.

[figure(s) omitted; refer to PDF]

The maximum singular value curve (Figure 23) of the corresponding signal is calculated from the atmosphere turbulence response of the four-channel. Figures 24–26 show the system identification results of the different analysis frequency bands, mainly including those at 3.3 Hz, 7 Hz, 17 Hz, and 32 Hz.

[figure(s) omitted; refer to PDF]

The system identification results illustrate the proposed method from the perspective of frequency domain data fitting. Table 3 analyzes the modal parameter calculation results of the frequency-spatial domain decomposition (FSDD) [14], SSI [15], and the proposed method.

Table 3

System identification results of the physical test.

	FSDD		SSI		Proposed method
	Frequency (Hz)	Damping ratio	Frequency (Hz)	Damping ratio	Frequency (Hz)	Damping ratio
Horizontal tail signal	10.351	4.461	10.573	6.871	11.396	9.973
	25.186	4.947	25.836	4.458	25.820	6.594
Wing signal	3.228	7.674	3.424	10.485	3.337	17.92
	6.521	3.861	6.832	6.083	6.774	4.754
	17.973	4.482	17.356	6.832	17.972	13.881
	32.217	4.634	32.638	5.489	32.622	5.764

The FSDD method is an improved frequency domain decomposition that estimates the modal parameters of random excitation signals through the frequency domain identification of the SDOF system. This approach is applied to analyze vibration responses in civil engineering [16] and system identification and modal extraction from response data [17]. The FSDD method estimates parameters in the frequency domain for a SDOF system and must decompose the maximum singular value by taking the frequency of interest as a parameter to separate the SDOF system. The SSI technique is a modal parameter estimation method based on state space modeling. For the output-only system identification problem, modal parameter estimation is performed based on the relationship between the system matrix and the Hankle matrix of the response signal. In studies involving the SSI method, it is usually necessary to estimate the modal parameters of multiple orders and further extract the stable modal through the steady-state diagram method. Because the modal parameter estimation of multiple orders is required, the algorithm execution time of the SSI is too long.

The modal parameters of the flutter test have no true values and are thus analyzed using FSDD, SSI, and the proposed method. The analysis indicates that the estimated results of the frequency of interest are essentially the same. Because the damping parameter is relatively small, there are cases where the error of individual estimated parameters is large. Table 4 compares the results of the three different methods based on the modal (17 Hz and 33 Hz) of the left and right symmetrical positions of the wing using the proposed method. The results show that the frequency estimates of the left and right symmetrical positions are generally consistent in terms of the data acquisition error. For the damping parameters, there are errors in the estimation results of the response signals at different locations.

Table 4

Modal parameter estimation results of the wing symmetrical position.

All wing signal		Left wing signal		Right wing signal
Frequency (Hz)	Damping ratio	Frequency (Hz)	Damping ratio	Frequency (Hz)	Damping ratio
17.972	13.881	17.374	10.581	18.918	12.114
32.622	5.764	33.098	5.340	33.269	5.226

The results show that the model parameters estimated using the proposed method are similar with those obtained using FSDD and SSI. The modal parameter estimation method for the symmetrical position of the wing also demonstrates the reliability of the proposed method. The maximum singular value curve based on the rational fraction polynomial can be directly used in the analysis frequency interval. The system identification is performed at the frequency band, and the efficiency of the modal parameter estimation is further improved. The frequency domain fitting error caused by the spectral pin density of the autocorrelation PSD function of the SDOF system is further avoided.

4.2. Flutter Boundary Prediction

The 3.3 Hz mode of the wing acceleration signal is taken as an example to illustrate the engineering applicability of the proposed method for flutter boundary prediction. The estimation of the modal parameters (frequency, damping ratio) in the boundary prediction is based on the turbulence response of the aircraft at a stable level flight under multiple speed steps. The modal parameter estimation of the response signal of each speed step is carried out using the proposed method, and the final frequency varies with the velocity shown in Figure 26. The modal parameter identification of the damping ratio result is carried out under the condition that the frequency identification is accurate. The boundary prediction result of the final damping ratio is shown in Figure 27, in which the flutter boundary velocity is close to 0.9734 Ma and thus consistent with the flutter phenomenon.

[figure(s) omitted; refer to PDF]

5. Conclusion

This paper analyzes the turbulence response characteristics of multichannels based on the randomness of atmosphere turbulence excitation and frequency domain orthogonality of the signal. The power spectral density function matrix of the turbulence response of multichannels is calculated, and singular value decomposition is performed for the power spectrum density function matrix at each spectral pin using the orthogonality of the frequency domain. The maximum singular value curve of the entire frequency band of the system can thus be calculated, which represents the auto power spectral density function of the system. The system identification for the maximum singular value curve using rational fractional polynomials can determine the frequency response function of the system. Further analysis can be performed to obtain the modal parameters of the system. A theoretical verification of the method is carried out using the simulation signal, and the flutter test data of a certain flight type are used for engineering verification, which verify the theoretical correctness and engineering applicability of the proposed method.

Acknowledgments

This work is supported by the Innovation Foundation for Doctor Dissertation of Northwestern Polytechnical University (Grant no. CX2021077) and the Fundamental Research Funds for the Central Universities (Grant no. 31020190MS702).

References

[1] O. Y. Al-Thehabey, "Dynamic instability analysis of aeroelastic systems with application to aircraft wings," Physics of Fluids, vol. 33 no. 9, article 094110,DOI: 10.1063/5.0063901, 2021.

[2] J. F. Kiviaho, K. E. Jacobson, G. J. Kennedy, "Flutter boundary identification from time-domain simulations using the matrix pencil method," AIAA Journal, vol. 57 no. 8, pp. 3639-3645, DOI: 10.2514/1.J058072, 2019.

[3] N. Pandiya, W. Desmet, "Direct estimation of residues from rational-fraction polynomials as a single- step modal identification approach," Journal of Sound and Vibration, vol. 517, article 116530,DOI: 10.1016/j.jsv.2021.116530, 2022.

[4] T. Wei, "Flutter modal identification with frequency-based least squares estimator using vector orthogonal polynomials," 2008 27th Chinese control conference, pp. 265-269, .

[5] Q. Wen, X. G. Hua, Z. Q. Chen, H. W. Niu, X. Y. Wang, "AMD-based random decrement technique for modal identification of structures with close modes," Journal of Aerospace Engineering, vol. 31 no. 5, article 04018057,DOI: 10.1061/(ASCE)AS.1943-5525.0000882, 2018.

[6] S. Oymak, N. Ozay, "Revisiting Ho-Kalman based system identification: robustness and finite-sample analysis," IEEE Transactions on Automatic Control, 2022.

[7] M. Scionti, J. Lanslots, I. Goethals, A. Vecchio, H. Van der Auweraer, B. Peeters, B. De Moor, "Tools to improve detection of structural changes from in-flight flutter data," Proceedings of the Eighth International Conference on Recent Advances in Structural Dynamics, .

[8] R. Brincker, L. Zhang, P. Andersen, "Modal identification of output-only systems using frequency domain decomposition," Smart Materials and Structures, vol. 10 no. 3, pp. 441-445, DOI: 10.1088/0964-1726/10/3/303, 2001.

[9] J. Zeng, S. L. Kukreja, "Flutter prediction for flight/wind-tunnel flutter test under atmospheric turbulence excitation," Journal of Aircraft, vol. 50 no. 6, pp. 1696-1709, DOI: 10.2514/1.C031710, 2013.

[10] Y. Govers, H. Mai, J. Arnold, J. K. Dillinger, A. K. Pereira, E. K. Takara, M. S. Correa, O. A. Mello, R. F. Marques, E. G. Geurts, R. J. Creemers, "Wind tunnel flutter testing on a highly flexible wing for aeroelastic validation in the transonic regime within the HMAE1 project," IFASD 2019-International Forum on Aeroelasticity and Structural Dynamics, .

[11] G. Jelicic, J. Schwochow, Y. Govers, J. Sinske, R. Buchbach, J. Springer, "Online monitoring of aircraft modal parameters during flight test based on permanent output-only modal analysis," 58th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, .

[12] J. S. Bendat, A. G. Piersol, Random Data: Analysis and Measurement Procedures, [M.S. thesis], 2010.

[13] M. H. Richardson, D. L. Formenti, "Parameter estimation from frequency response measurements using rational fraction polynomials," Proceedings of the 1st international modal analysis conference, pp. 167-186, .

[14] L. Zhang, T. Wang, Y. Tamura, "A frequency-spatial domain decomposition (FSDD) method for operational modal analysis," Mechanical Systems and Signal Processing, vol. 24 no. 5, pp. 1227-1239, DOI: 10.1016/j.ymssp.2009.10.024, 2010.

[15] V. Boonyapinyo, T. Janesupasaeree, "Data-driven stochastic subspace identification of flutter derivatives of bridge decks," Journal of Wind Engineering and Industrial Aerodynamics, vol. 98 no. 12, pp. 784-799, DOI: 10.1016/j.jweia.2010.07.003, 2010.

[16] S. T. K. Lin, Y. Lu, M. M. Alamdari, N. L. D. Khoa, "Field test investigations for condition monitoring of a concrete culvert bridge using vibration responses," Structural Control and Health Monitoring, vol. 27 no. 10,DOI: 10.1002/stc.2614, 2020.

[17] P. C. Schulze, P. Thompson, B. P. Danowsky, D. C. Lee, C. Harris, "System identification and modal extraction from response data," AIAA Atmospheric Flight Mechanics (AFM) Conference, .