# Integration of Freeplay-Induced Limit Cycles Based On a State Space Iterating Scheme
## 4 Numerical Model and Results of the Henon-RK45 Method  
In this section, a plunge-pitch wing section is introduced, as shown in Figure 9. Time integration with the Henon-RK45 method is used to obtain the nonlinear LCO behaviors when symmetrical and non-preloaded freeplay nonlinearity is added to the pitch DoF of the wing section. The amount of freeplay is \(2\delta=0.1^{\circ}\). The results in this section will serve as the benchmark for the next section, where the SSI scheme is applied to the same numerical model. The parameters of the wing section are listed in Table A1 in Appendix A.  
Figure 9: The numerical model: a plunge–pitch wing section, where freeplay nonlinearity is introduced in the pitch degree of freedom (DoF). An unsteady aerodynamic model is established in the frequency domain following [PERSON]'s method [27] and can be expressed as \(\mathbf{f}_{a}(k)=\frac{1}{2}\rho UU^{2}\mathbf{Q}\mathbf{q}(k)\), where \(\rho=1.225\,\mathrm{kg/m^{3}}\) is the air density, \(k\) is the reduced frequency, \(U\) is the flow velocity (m/s), \(\mathbf{Q}\) is the aerodynamic influence coefficient (AIC) matrix, and \(\mathbf{q}=\left[h,\ \alpha\right]^{T}\) are the generalized modal coordinates, where \(h\) and \(\alpha\) denote the plunge (m) and pitch angle (rad), respectively. Then, the rational function approximation (RFA) introduced by [PERSON] [28] and developed by [PERSON] [29] and the minimal state (MS) method proposed by [PERSON] and [PERSON] [30] are used to approximate the aerodynamic loads in the time domain, \(\mathbf{\tilde{f}}_{a}(t)\). The equations of motion for the wing section model are as follows:  
\[ \begin{split}(\mathbf{M}_{\mathrm{S}}\mathbf{\ddot{q}}+\mathbf{ D}_{\mathrm{S}}\mathbf{\dot{q}}+\mathbf{K}_{\mathrm{S1}}\mathbf{q})/L= \mathbf{\tilde{f}}_{a}(t)&=\frac{1}{2}\rho U^{2}\mathbf{A}_{ 0}\mathbf{q}+\frac{1}{2}\rho bU\mathbf{A}_{1}\mathbf{\dot{q}}+\frac{1}{2} \rho b^{2}\mathbf{A}_{2}\mathbf{\ddot{q}}+\frac{1}{2}\rho U^{2}\mathbf{A}_{ D}\mathbf{r}\\ \mathbf{\dot{r}}&=\mathbf{A}_{E}\mathbf{\dot{q}}+U/b \mathbf{R}\end{split}  \tag{13}\]  
where \(\mathbf{M}_{\mathrm{S}}\), \(\mathbf{D}_{\mathrm{S}}\), and \(\mathbf{K}_{\mathrm{S1}}\) represent the structural mass, damping, and stiffness, respectively; \(\mathbf{A}_{0}\), \(\mathbf{A}_{1}\), \(\mathbf{A}_{2}\), \(\mathbf{A}_{D}\), and \(\mathbf{A}_{E}\) are matrices of the aerodynamic coefficients; and \(\mathbf{R}\) is a diagonal matrix that contains two lag terms, \(r_{1}=-0.08\) and \(r_{2}=-0.60\). The accuracy of \(\mathbf{\tilde{f}}_{a}(t)\) approximating \(\mathbf{f}_{a}(k)\) in Equation (13) is examined by converting \(\mathbf{\tilde{f}}_{a}(t)\) in the Laplace domain, where \(L(\mathbf{\tilde{f}}_{a}(t))=\)\(\frac{1}{2}\rho U^{2}\mathbf{\dot{Q}}(s)\mathbf{q}(s)\). Next, let the Laplace operator be \(s=ik\), where \(i\) is an imaginary unit. Figure 10 shows a comparison between \(\mathbf{\ddot{Q}}(ik)\) and \(\mathbf{Q}(k)\) with 18 reduced frequencies: \(k=0.005\), \(0.01\), \(0.02\), \(0.03\), \(0.04\), \(0.05\), \(0.06\), \(0.08\), \(0.10\), \(0.12\), \(0.15\), \(0.20\), \(0.25\), \(0.30\), \(0.35\), \(0.40\), \(0.45\), and \(0.50\). This result proves that the aerodynamic model in Equation (13) is qualified for further calculations. When freeplay nonlinearity is added to the pitch DoF, the linear structurally restoring forces \(\mathbf{K}_{\mathrm{S1}}\mathbf{q}\) will be replaced by \(\mathbf{K}_{\mathrm{S0}}\mathbf{q}\).