## 2 Aeroelastic Model and Analysis Method
### Analysis Methods for Aeroelastic Model with Freeplay
#### 2.3.2 Numerical Integration Approach  
According to Equations (5) and (11), the state-space equations of the nonlinear 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}\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}+\begin{bmatrix}\mathbf{0}\\ -\mathbf{M}^{-1}\mathbf{F}_{\mathrm{NL}}\end{bmatrix} \tag{14}\]where  
\[\mathbf{F}_{\text{NL}}=\begin{bmatrix}M_{\text{B}}\\ M_{\text{T}}\end{bmatrix}\]  
In this study, the Henon-RK45 method [13; 14] is used to integrate Equation (14). For simplicity, Equations (14) can be expressed as \(\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})\). The fourth-fifth-order Runge-Kutta scheme (RK45) is used for numerical integration until the monitored variable \(\mathbf{x}^{(k)}\) passes through the upper or lower boundary of freeplay. Assume that the value of the monitored variable \(\mathbf{x}^{(k)}\) is \(\mathbf{x}_{c}^{(k)}\) at the crossing point (\(\mathbf{x}_{c}^{(k)}\) is \(\delta\) or \(-\delta\)), the last time instance before crossing is \(t_{0}\), and the corresponding states are \(\mathbf{x}_{0}\). Then, the Henon method is used to make the monitored variable \(x_{k}\) reach the discontinuity point \(\mathbf{x}_{c}^{(k)}\) precisely via the following steps:  
1. change the integral variable from time \(t\) to the monitored variable \(\mathbf{x}^{(k)}\), and assume a new vector of the states \(\mathbf{y}=[\mathbf{x}^{T},t]^{T}\); that is, \[\frac{\mathbf{d}\mathbf{y}}{\text{d}\mathbf{x}^{(k)}}=\frac{\mathbf{d} \mathbf{y}}{\text{d}t}\cdot\frac{\text{d}t}{\text{d}\mathbf{x}^{(k)}}=\begin{bmatrix} \mathbf{f}(\mathbf{x})/\mathbf{f}^{(k)}(\mathbf{x})\\ 1/\mathbf{f}^{(k)}(\mathbf{x})\end{bmatrix}\]
2. integrate the above equation using the RK45 method from \(\mathbf{y}_{0}=[\mathbf{x}_{0}^{T},t_{0}]^{T}\) to the crossing point \(\mathbf{y}_{c}=[\mathbf{x}_{c}^{T},t_{c}]^{T}\) as the integral variable \(\mathbf{x}^{(k)}\) is moved from \(\mathbf{x}_{0}^{(k)}\) to \(\mathbf{x}_{c}^{(k)}\);
3. restore the integral variable to the time \(t\), and integrate Equations (14) from \(\mathbf{x}(t_{c})=\mathbf{x}_{c}\) using the RK45 method until the next crossing point. Then repeat steps (1) to (3) until the end.  
For systems with a single freeplay nonlinearity, the Henon-RK45 method can effectively identify the discontinuities and achieve the purpose of reducing the error. For systems with multiple freeplay nonlinearities, there may be two or more integral variables that simultaneously cross the upper or lower boundary of the freeplay in the same time step. In this case, it is necessary to determine the order in which the relevant integral variables arrive at the boundary, and then the integral variable and time step can be controlled so that the state of the system can reach these boundaries successively.  
Regarding the integration time step, it is recommended to set it around \(1/(20\omega_{\text{max}})\), where \(\omega_{\text{max}}\) is the maximum modal frequency of the model modes considered.