## 2 Aeroelastic Model and Analysis Method
### Analysis Methods for Aeroelastic Model with Freeplay  
#### 2.3.1 Describing Function Method  
Common linear flutter solution methods, such as the V-g method, p-k method, and root-locus method, have high computational efficiency. However, when freeplay nonlinearity is involved, these methods are no longer applicable. In this instance, an equivalent linearization method, the DF method, has solved this problem. Taking into account only the first-order harmonic component of the output of the nonlinear element under the harmonic input, the relationship between the equivalent stiffness and the LCO amplitude is established as follows [5]:  
\[\frac{\bar{K}}{K}=\begin{cases}0&A\leq\delta\\ 1-\frac{2t_{1}}{\pi}-\frac{\sin 2t_{1}}{\pi}&A>\delta\end{cases} \tag{12}\]  
where \(t_{1}=\arcsin(\delta/A)\), in which \(\delta\) is the unilateral freeplay size and \(A\) is the LCO amplitude in the nonlinear element. \(A\leq\delta\) indicates that the system has been vibrating inside the freeplay in the nonlinear element, that is, the elastic restoring moment and the equivalent stiffness are 0. When \(A>\delta\), the equivalent stiffness increases from 0 to 1 with increasing oscillatory amplitude. Here, we define the amplitude ratio \(A/\delta\) as the nondimensional LCO amplitude.  
According to Equations (3) and (11), the state-space equations of the linear system can be written in the following form:  
\[\begin{bmatrix}\dot{\mathbf{q}}\\ \dot{\mathbf{q}}\end{bmatrix}=\begin{bmatrix}\mathbf{0}&\mathbf{I}\\ \mathbf{M}^{-1}\Big{(}\frac{1}{2}\rho V^{2}\mathbf{A}_{0}-\mathbf{K}\Big{)}& \mathbf{M}^{-1}\Big{(}\frac{1}{2}\rho V\mathbf{A}_{1}\Big{)}\end{bmatrix} \begin{bmatrix}\mathbf{q}\\ \dot{\mathbf{q}}\end{bmatrix} \tag{13}\]  
where  
\[\mathbf{M}=\begin{bmatrix}I_{xx}&-I_{xy}\\ -I_{xy}&I_{yy}\end{bmatrix},\mathbf{K}=\begin{bmatrix}K_{\mathrm{B}}\\ &K_{\mathrm{T}}\end{bmatrix},\mathbf{q}=\begin{bmatrix}q_{\mathrm{B}}\\ q_{\mathrm{T}}\end{bmatrix}\]  
Here, we define the matrix \(\mathbf{A}\) as follows:  
\[\mathbf{A}=\begin{bmatrix}\mathbf{0}&\mathbf{I}\\ \mathbf{M}^{-1}\Big{(}\frac{1}{2}\rho V^{2}\mathbf{A}_{0}-\mathbf{K}\Big{)}& \mathbf{M}^{-1}\Big{(}\frac{1}{2}\rho V\mathbf{A}_{1}\Big{)}\end{bmatrix}\]  
For systems with a single freeplay nonlinearity, when the nondimensional LCO amplitude in the nonlinear element is given, an equivalent stiffness \(\bar{K}\) can be obtained from Equation (12). By replacing \(K_{\mathrm{B}}\) or \(K_{\mathrm{T}}\) in Equation (13), a linear flutter velocity can be calculated through the root-locus method. The specific process is as follows:  
1. Assume a nondimensional LCO amplitude \(A\) in the nonlinear element and obtain the equivalent stiffness \(\bar{K}\);
2. Replace \(K_{\mathrm{B}}\) or \(K_{\mathrm{T}}\) with \(\bar{K}\) in Equation (13) and obtain a new matrix \(\mathbf{A}\);
3. Calculate the eigenvalues of \(\mathbf{A}\) with the variation of the free-stream velocity \(V\);
4. As \(V\) increases, find the eigenvalue first crossing the imaginary axis from negative to positive, and record the free-stream velocity \(V\), the crossing eigenvalue and corresponding eigenvector at this time;
5. Change the assumed nondimensional LCO amplitude and repeat steps (1) to (5) until the end.  
Finally, we can obtain the relationship between the free-stream velocity and the nondimensional LCO amplitude; this reflects the nonlinear aeroelastic characteristics of the system. The frequency and phase information of the LCO will be described in more detail in Section 2.3.3. For systems with multiple freeplay nonlinearities, a discrete scanning method will be introduced in Section 3.