## 3 Numerical Method Verification
### Simulation of AGARD445.6 Wing Flutter  
In the previous section, the dynamic model was verified by a single pendulum example. In order to verify the correctness of the CFD-based rigid-flexible coupling research method, the flutter boundary of the AGARD445.6 wing was predicted and compared with the experimental results. The AGARD445.6 wing has been frequently used to verify the reliability of the proposed method for aeroelastic problems. This model was used by Langley Research Center to study the wing flutter characteristics in transonic conditions based on public wind tunnel test data [28]. The spanwise section of the wing is based on NACA65A004 airfoil, with an aspect ratio of 1.644, root tip ratio of 0.659, and 1/4 chord sweep angle of 45\({}^{\circ}\). The wing model geometry is depicted in Figure 7.  
The modal characteristics of the AGARD 445.6 wing were investigated through structural modal analysis. The vibration behavior is primarily characterized by its first four dominant mode shapes: first bending mode, first torsion mode, second bending mode, and second torsion mode. To visualize the modal analysis results of the wing's structural model, these fundamental modes were interpolated onto the aerodynamic surface mesh and are comprehensively illustrated in Figure 8.  
The Mach numbers in computation are 0.499, 0.678, 0.901, 0.96, 1.072, and 1.141. The angle of attack of the incoming flow is 0\({}^{\circ}\), and the critical flutter velocity at each Mach number is determined by changing the incoming flow pressure. The dimensionless flutter velocity is expressed as follows:  
\[V*=\frac{V}{b_{s}\omega_{\alpha}\sqrt{\#}} \tag{29}\]  
where \(V\) represents the incoming flow velocity, \(b_{s}\) is the half chord length of the wing root, and \(\omega_{\alpha}\) denotes the first-order uncoupled torsional frequency. Prior to the computation, grid independence verification was conducted by employing three meshes with different sizes, as illustrated in Figure 9.  
Numerical simulations were performed at Mach 0.678 to calculate the drag coefficient of the AGARD445.6 wing under converged conditions for each grid configuration. The results presented in Figure 10 demonstrate that the drag coefficient obtained with the medium-density mesh showed minimal discrepancies compared to that of the refined mesh. To optimize computational efficiency while maintaining accuracy, the medium-density mesh was selected for subsequent aeroelastic simulations. Figure 7: The AGARD445.6 wing model geometry. Figure 8: The first four mode shapes of the AGARD445.6 wing. Figure 10: Drag coefficient results under different meshes. Figure 9: Surface meshes with varying sizes. Before commencing the calculation, an appropriate initial velocity must be applied to the first-order mode to serve as an excitation input, which is generally taken as \(\dot{q}=2\pi f_{1}q_{1}\), where \(f_{1}\) is the first-order modal frequency and \(q_{1}\) is a small generalized displacement, which can be taken as 0.001. The computational time step is 0.0002. The time-dependent variation in generalized displacement at Mach 0.678 under different dimensionless velocities is illustrated in Figure 11. As the figure shows, the system exhibits three distinct dynamic regimes: when dimensionless velocity \(V*=0.41\), the generalized displacement converges; when dimensionless velocity \(V*=0.4175\), the generalized displacement oscillates with equal amplitude; when the dimensionless velocity \(V*=0.42\), the generalized displacement diverges. The numerically predicted dimensionless flutter critical velocity of 0.4175 shows remarkable agreement with experimental measurements (0.4174), with a relative deviation of merely 0.024%. For each Mach number, the AGARD 445.6 wing corresponds to a specific flutter speed index. Figure 12 presents a comprehensive comparison of flutter speeds across six Mach numbers, integrating findings from the current study with existing reference data.