## 3 Numerical Validation
### FSI with Free Motion  
In this section, flutter prediction and analysis are conducted on the NACA64A010 airfoil and AGARD 445.6 wing. Figure 16 presents their configurations, respectively. For the sake of problem description, exemplified by the typical section of the wing model, AGARG 445.6 wing in 3D problem is viable to analogy. In light of Figure 16a, the model contains two degrees of freedom, referred to as pitching \(\alpha\) and plunging \(h\). The motions are governed by  
\[ \begin{split} m\ddot{h}+S_{\alpha}\ddot{\alpha}+K_{h}h& =-L\\ S_{\alpha}\ddot{h}+I_{\alpha}\ddot{\alpha}+K_{\alpha}\alpha& =M_{EA}\end{split}  \tag{54}\]  
where \(L\) and \(M_{EA}\) denote the lift and moment about the elastic axis, respectively. \(K_{h}\) and \(K_{\alpha}\) are blending and torsional spring stiffness. m represents the mass per unit span, \(S_{\alpha}\) denotes the static moment around the elastic axis, and \(I_{\alpha}\) indicates the rotational moment of inertia. With dimensionless time reduced by the uncoupled natural frequency of the torsional spring, \(\tau=\omega_{\alpha}t\), the equation of motion is rewritten as  
\[\mathbf{M}\ddot{\mathbf{q}}+\mathbf{K}\mathbf{q}=\mathbf{F} \tag{55}\]  
where  
\[\mathbf{M}= \begin{bmatrix}1&x_{\alpha}\\ x_{\alpha}&r_{\alpha}^{2}\end{bmatrix} , \tag{56}\]  
\[\mathbf{K}= \begin{bmatrix}\left(\frac{\omega_{\alpha}}{\omega_{\alpha}} \right)^{2}&0\\ 0&r_{\alpha}^{2}\end{bmatrix}  \tag{57}\]  
are the dimensionless mass and stiffness matrices,  
Figure 16: Schematic of (**a**) the typical section wing model and (**b**) AGARD 445.6 wing. \[\mathbf{F}=\frac{1}{\pi\mu k_{c}^{2}} \begin{bmatrix}-C_{l}\\ 2C_{m}\end{bmatrix} =\frac{U^{*2}}{\pi\mu} \begin{bmatrix}-C_{l}\\ 2C_{m}\end{bmatrix} , \tag{58}\] \[\mathbf{q}= \begin{bmatrix}\frac{h}{\theta}\\ \alpha\end{bmatrix}  \tag{59}\]  
states the combination of loads and displacements. \(C_{l}\) and \(C_{m}\) are the coefficients of lift and moments about the elastic axis. \(x_{\alpha}=S_{\alpha}/b\) is static unbalance, \(\omega_{h}=\sqrt{K_{h}/m}\) is the uncoupled natural frequency in plunge, \(\omega_{\alpha}=\sqrt{K_{\alpha}/I_{\alpha}}\) is the uncoupled natural frequency in pitch. \(r_{\alpha}^{2}=I_{\alpha}/\left(m\theta^{2}\right)\) is the squared radius of gyration, \(U^{*}\) is the reduced velocity defined as \(U_{\infty}/\left(\omega_{\alpha}b\right)\), \(k_{c}=\omega_{\alpha}b/U_{\infty}\) is the reduced frequency, and \(\mu=m/\left(\rho_{\infty}v\right)\) is the mass ratio. \(\rho_{\infty}\) is endowed with the density of flow and \(v\) is the volume coefficient. As a matter of fact, a dimensionless parameter, flutter velocity \(V_{f}^{*}\), is put forth to evaluate the critical value. When local velocity exceeds that, a small disturbance is going to amplify. \[V_{f}^{*}=\frac{U_{\infty}}{b\omega_{\alpha}\sqrt{\mu}} \tag{60}\]  
The function of \(V_{f}^{*}\) depends on a Mach number, so there exists a flutter boundary differing the situations into stable and unstable conditions. From a practical point of view, good prediction about flutter boundary plays an important role in the transonic region. The indicator \(V_{f}^{*}\) measures the flutter point in quantity. When \(V_{f}^{*}\) is lower than the critical value, the oscillation is damped. As \(V_{f}^{*}\) increases, the neutrally stable point is reached. When \(V_{f}^{*}\) exceeds the critical value, the amplitude of oscillation diverges, even inducing collapse. To capture the flutter boundary, a series of simulations are performed by varying the \(V_{f}^{*}\). In a standard rule, various runs are needed until a zero damping response is yielded. However, the process of operation in this work is simplified under the guidance of the literature [56]. Owing to their accessible flutter boundary, for each Mach number, only several \(V_{f}^{*}\) spanning beside the sensible and sensitive regions are simulated, with an increment \(\Delta V_{f}^{*}=0.1\).