## 2 Numerical and Experimental Methods
### Vortex Particle Simulation (VPM)
\begin{table}
\begin{tabular}{l c c c c} \hline \hline Reynolds number, \(Re_{C}\) & 4000 & & 60,000 \\ \hline Time increment, \(\Delta t\) (sec) & 0.005 & & 0.002 \\ \hline Particle spacing at the center of the airfoil, \(h\) (\(\times 10^{-3}\) m) & 1.6 & & 0.25 \\ \hline Airfoils & NACA64 & DU40 & NACA64 & DU40 \\ \hline Surface panel length, \(\Delta S_{B}\) (\(\times 10^{-3}\) m) & 0.7–3.4 & 0.6–1.8 & 0.4–1.7 & 0.32–0.9 \\ \hline Azimuthal division of the redistribution grid, \(N_{m}\) & 8024 & & 49217 \\ \hline Control volume for integral force calculation, Equation (12) & & & \\ \hline Control volume limits in x-dir, \([x_{\text{min}},x_{\text{max}}]\) (m) & & \([-0.8,1.5]\) & \\ \hline Control volume limits in y-dir, \([y_{\text{min}},y_{\text{max}}]\) (m) & & \([-1.2,1.2]\) & \\ \hline Eulerian grid for pressure Poisson solution, Equation (13) & & & \\ \hline Domain limits in x-dir, \([x_{\text{min}},x_{\text{max}}]\) (m) & \([-1.5,5]\) & \([-0.8,1.3]\) & \\ \hline Domain limits in y-dir, \([y_{\text{min}},y_{\text{max}}]\) (m) & \([-3,3]\) & \([-0.8,0.8]\) & \\ \hline Grid spacing, \(\Delta S_{G}\) (\(\times 10^{-3}\) m) & 1.6 & 0.26 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Numerical setup in the vortex particle method (VPM). Figure 3: Time evolutions of flow quantities: (**a**) Numbers of particles used in simulations. (**b**) Lift force. (**c**) Drag force. The integral forces, \(\mathbf{F}=F_{D}\mathbf{e}_{1}+F_{L}\mathbf{e}_{2}\), are calculated using the integral momentum equations derived by [28] that are convenient for the current VPM:  
\[\mathbf{F}=-\frac{d}{dt}\int_{V}\rho\,\mathbf{x}\times\mathbf{\omega}\ dV+\oint_{S} \rho\,\mathbf{n}\cdot\mathbf{\gamma}\ dS \tag{12}\]  
where the lift, \(F_{L}\), and drag, \(F_{D}\) are defined in Figure 1. Additionally, \(\mathbf{\gamma}=0.5\mu^{2}\mathbf{I}-\mathbf{u}\mathbf{u}-\mathbf{u}(\mathbf{x} \times\mathbf{\omega})+\mathbf{\omega}(\mathbf{x}\times\mathbf{u})\) is the tensor that needs to be evaluated on the surface, \(S\), of the control volume, \(V\). The viscous components of the force are neglected in the original formulation due to the small portion of the overall force (generally of a few percent, see, e.g., [6]). Furthermore, since the airfoil body is assumed to be stationary, the control volume and its surface remain fixed. The setup of the control volume for the force calculation is shown in Table 1. The resulting lift and drag force histories are shown in Figure 3b,c, respectively, with legends consistent with those of Figure 3a. After the impulsive start of the airfoils, traces of force development can be observed before the steady-state fluctuations occur. The steady-state portions of the force histories are used to obtain the representative force statistics. For simulations of \(Re_{C}=4000\), the steady-state forces generally occur for \(5\leq Ut/C\leq 50\), leading to time-averaging durations of \(UT_{\text{avg}}/C=45\). For simulations of \(Re_{C}=60\),000, the steady-state forces occur for \(10\leq Ut/C\leq 15\), giving the time-averaging durations of \(UT_{\text{avg}}/C=5\). Moreover, instants of pressure fields are calculated by solving the Poisson equations of pressure, for example, as follows:  
\[\
abla^{2}p=\rho\
abla\cdot\bigg{(}\mathbf{u}\times\mathbf{\omega}-\
abla\bigg{(} \frac{u^{2}}{2}\bigg{)}-v\
abla\times\mathbf{\omega}\bigg{)}. \tag{13}\]  
where the velocity gradient sources, such as \(\mathbf{\omega}\) and \(\
abla\big{(}u^{2}/2\big{)}\), are'measured' on the Eulerian grids. The solution is obtained on a rectangular domain enclosing the airfoil. Neumann boundary conditions are used on the outlet border while Dirichlet conditions are imposed on the rest of the boundaries. The Dirichlet conditions are obtained using Bernoulli's equation.