## 2 Methodology and Validation  
This section presents the development and validation of the numerical method, which is established to handle blade rotation problems. Under the framework of an in-house code based on the finite volume method, the third-order scheme of MUSCL-Roe [16] is implemented first for the spatial discretization of the convective fluxes, and then the viscous fluxes are discretized with the second-order central differences. The dual-time approach based on the fully implicit LU-SGS method [17] is implemented as the time-marching method. The effect of turbulence is estimated using the Spalart-Allmaras one-equation model [18]. The reliability of the code has been verified via a series of works [19; 20]. The feasibility of the numerical method for the numerical investigation is verified through the case of a two-bladed Caradonna-Tung rotor under the forward flight condition [21]. The purpose of the test is to obtain the pressure distribution and aerodynamic load on the rotor surface. The radius of the rotor is 1.143 m and the aspect ratio of the blade is 6. The cross-sectional shape is the NACA0012 airfoil. Pressure sensors are installed on two blades in the same radial positions. The ICEM-Hexa is applied in the generation of a multi-block structured grid. The O-topology is applied to describe the rotation flow field. The grid topology of the rotor is shown in Figure 1. For the numerical simulation issues with relative motions, the dynamic patched grid is used [22; 23], which has the primary advantage of handling relative motion. Based on this grid generation method, the grid of the helicopter and main rotor can be generated independently without considering the topology relation among adjacent blocks. Fluxes are interpolated through the patched surfaces, thus enabling the coupling of the flow information of adjacent subdomains. For the purpose of improving the prediction accuracy of the rotor wake vortex and flow field, the grid in the region below the rotor is refined. The first cell height near the wall is set at 1 \(\times\) 10\({}^{-6}\) R to ensure y\(+\)\(\approx\) 1, and the complete grid consists of 10,489,523 nodes. Figure 2 demonstrates the spatial grid of the blades and the grid used for the boundary layer profile. The numerical simulation is performed for the forward flight case, corresponding to the conditions shown in Table 1. In the current study, the step size of the per physical time step is set at equivalent to \(\Delta\Psi=1.0\) deg. The numerical simulation is carried out for 50 full revolutions first to provide a convergent flow field solution, then for 10 full revolutions to obtain a reliable solution, which is consistent with previous works [20]. The total numerical effort in CPU time is 200 h and the number of processors is 64. Figure 3 shows the comparison between the calculated chordwise pressure distribution and the test result at the radial position \(r/R=0.89\) at several azimuth angles in one rotor revolution. It can be seen that the numerical results are in fairly good agreement with the experimental results at most azimuth angles. The pressure peak at the leading edge of the \begin{table}
\begin{tabular}{c c} \hline Main rotor advanced ratio & 0.3 \\ Freestream Mach number & 0.188 \\ Main rotor tip Mach number & 0.628 \\ Reynolds number & \(3.66\times 10^{6}\) \\ \hline \end{tabular}
\end{table} 
Table 1: Computational conditions for the forward flight case. Figure 1: Topology of the computational domain. Figure 2: Spatial grid of the blades and the grid used for the boundary layer profile. upper surface of the blade and the subsequent pressure recovery process are predicted well. However, noticeable discrepancies can be seen at \(\Psi=135^{\circ}\), mainly because the strength of the shock wave near the radial station is notable at this particular azimuth angle. The isosurface of the Q-criterion is shown in Figure 4. The result shows that vortices induced by the blade tip extend helically backwards and interfere with each other in the forward flight case.