## 3 Numerical Simulation Method  
In order to generate the computational grids around the OA209 airfoil, the Poisson equations [28] are chosen as the governing equations under 2D conditions. By solving these equations, the C-topology grids with \(459\times 80\) points are generated. To satisfy the requirement of OFV simulation, the far-field boundary of airfoil grids is fixed at 25 times the airfoil chord, and the y\({}^{+}\) of the grids near the wall surface is smaller than 1.2. In this research, the background grid is composed of \(300\times 300\) points with a far-field boundary of 50 c (airfoil chord), and the inverse map method [29] is employed in this code to search donor element. As shown in Figure 2, the oscillating velocity is simulated by moving the airfoil back and forth periodically [30; 31]. Figure 1: Mach number distribution of rotor disk and angle of attack (AoA) variation with rotor azimuth. (**a**) Mach number distribution. (**b**) AoA and freestream velocity variations. In order to simulate the unsteady compressible flowfield around an airfoil, the integral form of the Navier-Stokes equations is employed in this work, i.e.,  
\[\frac{\partial}{\partial t}\iiint\limits_{\Omega}\mathbf{W}d\Omega+\iint\limits_ {\partial\Omega}\left(\mathbf{F}_{c}-\mathbf{F}_{v}\right)d\mathbf{S}=0 \tag{4}\]  
where \(\mathbf{W}\) represents the conserved variables, \(\mathbf{F}_{c}\) and \(\mathbf{F}_{v}\) denote convective fluxes and viscous fluxes, respectively, i.e.,  
\[\mathbf{W}=\left[ \begin{array}{c}\rho\\ \rho u\\ \rho v\\ \rho E\end{array} \right]\mathbf{F}_{c}=\left[ \begin{array}{c}\rho\mathbf{V} _{r}\\ \rho u\mathbf{V}_{r}+n_{x}p\\ \rho v\mathbf{V}_{r}+n_{y}p\\ \rho H\mathbf{V}_{r}+\mathbf{V}_{t}p\end{array} \right]\mathbf{F}_{v}=\left[ \begin{array}{c}0\\ n_{x}\tau_{xx}+n_{y}\tau_{xy}\\ n_{x}\tau_{yx}+n_{y}\tau_{yy}\\ n_{x}\Theta x_{x}+n_{y}\Theta y\end{array} \right] \tag{5}\]  
where \(\mathbf{V}_{r}=\mathbf{V}-\mathbf{V}_{t}\), \(\mathbf{V}\) is absolute velocity, \(\mathbf{V}_{t}\) is contravariant velocity. \(\tau_{ij}\) denotes the viscous stresses, and \(\Theta_{i}\) is the term describing the heat condition in the fluid. Meanwhile, the unsteady condition simulation is accomplished by employing the dual time-stepping approach. Then Equation (4) would be changed as  
\[\frac{d(\Omega^{n+1}\mathbf{W}^{*})}{d\tau}+\mathbf{R}^{*}(\mathbf{W}^{*})=0 \tag{6}\]  
where \(\mathbf{W}^{*}\) is the approximation for \(\mathbf{W}^{n+1}\), \(\tau\) represents a pseudo time variable, and the unsteady residual is defined as  
\[\mathbf{R}^{*}(\mathbf{W}^{*})=\frac{3\mathbf{W}^{*}-4\mathbf{W}^{n}+\mathbf{ W}^{n-1}}{2\Delta t}+\mathbf{R}(\mathbf{W}^{*}) \tag{7}\]  
To improve the computational efficiency, the implicit LU-SGS [32; 33] scheme is employed in this work. The two equation SST k - \(\omega\) turbulence model is employed to simulate the viscous stresses of the flowfield. In order to verify the accuracy of the present grid system and flowfield solver, a case of the NACA4412 airfoil is presented in this section. In this case, the pressure coefficient (Cp) simulated by the present CFD method compared with the test data [34] is shown in Figure 3a. In this case, three different grids (\(699\times 100\), \(459\times 80\), \(319\times 60\)) are used to analyze the grid sensitivity. By comparing with the test data, it illustrated that the Cp simulated by the CFD method correlates well with test data for a different grid. However, the convergence of a larger grid is lower than coarse grid, but it still converges at least five orders of magnitude. It is indicated that the grid used in this research could satisfy the calculation requirement. Meanwhile, it is illustrated in Figure 3b that the velocity profiles are also close to the test data. As a result, it is indicated that the grid used in this research is suitable. Figure 2: Schematic of moving-embedded grids.