## 2 Computational Model  
### Numerical Method  
The 2D, incompressible URANS equations are solved by ANSYS/FLUENT, for simulating the unsteady flow around a pitching airfoil. The governing equations are written as:  
\[\frac{\partial\rho}{\partial t}+\frac{\partial(\rho u_{i})}{\partial x_{i}}=0 \tag{1}\]  
\[\frac{\partial(\rho u_{i})}{\partial t}+\rho u_{j}\frac{\partial u_{i}}{ \partial x_{i}}=-\frac{\partial p}{\partial x_{i}}+\rho v\
abla^{2}u_{i}+F_{ P,i} \tag{2}\]  
where \(x_{i}\) stands for a coordinate direction, \(u_{i}\) denotes a component of fluid velocity vector, \(\rho\), \(p\), and \(v\) are the fluid density, pressure, and kinematic viscosity coefficient, respectively, and \(F_{P,i}\) refers to the momentum source term, accounting for the particle effect on the fluid phase.  
According to the conclusion by [PERSON] and [PERSON] [20], considering the transition from laminar to turbulent flow is crucial to accurately predict dynamic stall. Therefore, the transition SST model (also known as the \(\gamma\)-Re\({}_{\theta}\) model) is adopted for turbulence modeling, which couples the SST \(k\)-\(\omega\) eddy viscosity model [35] and two other transport equations [36].  
\[\frac{\partial k}{\partial t}+\frac{\partial}{\partial x_{j}}(u_{j}k)=\frac{ 1}{\rho}\frac{\partial}{\partial x_{j}}\left(\Gamma_{k}\frac{\partial k}{ \partial x_{j}}\right)+\frac{G_{k}-Y_{k}}{\rho} \tag{3}\]\[\frac{\partial\omega}{\partial t}+\frac{\partial}{\partial x_{j}}\big{(}u_{j} \omega\big{)}=\frac{1}{\rho}\frac{\partial}{\partial x_{j}}\Bigg{(}\Gamma_{ \omega}\frac{\partial\omega}{\partial x_{j}}\Bigg{)}+\frac{G_{\omega}-Y_{k}+D _{\omega}}{\rho} \tag{4}\]  
\[\frac{\partial\gamma}{\partial t}+\frac{\partial(u_{j}\gamma)}{\partial x_{j}} =\frac{P_{\gamma 1}-E_{\gamma 1}+P_{\gamma 2}-E_{\gamma 2}}{\rho}+\frac{1}{\rho} \frac{\partial}{\partial x_{j}}\Bigg{[}\bigg{(}\mu+\frac{\mu_{t}}{\sigma_{ \gamma}}\bigg{)}\frac{\partial\gamma}{\partial x_{j}}\Bigg{]} \tag{5}\]  
\[\frac{\partial\text{Re}_{\phi t}}{\partial t}+\frac{\partial\big{(}u_{j}\text{ Re}_{\phi t}\big{)}}{\partial x_{j}}=\frac{P_{\phi t}}{\rho}+\frac{1}{\rho} \frac{\partial}{\partial x_{j}}\Bigg{[}\sigma_{\theta t}(\mu+\mu_{t})\frac{ \partial\text{Re}_{\phi t}}{\partial x_{j}}\Bigg{]} \tag{6}\]  
The further description of the variables in these four equations can be found in the cited papers. This model has been proven to be reliable in predicting dynamic stall of airfoils [28; 31; 32; 37].  
For the simulation of unsteady incompressible flows, the PISO algorithm [38] is suitable, which has been proven by a series of research [28; 31; 32; 37]. In the spatial discretization, the least square cell-based and second-order are employed to compute the gradient and pressure equation. The other convection-diffusion equations are solved with a third-order MUSCL scheme. A second-order implicit scheme is used for time-dependent solution formulation.  
The discrete phase model (DPM) is employed to calculate the particle motion. A two-way coupled scheme is adopted, and particles are assumed to be non-interacting, non-deforming, and non-evaporating spheres. Then, the particle trajectories can be solved by:  
\[\frac{dv_{p}}{dt}=\beta\big{(}v_{a}-v_{p}\big{)}+g\bigg{(}\frac{\rho_{p}-\rho _{a}}{\rho_{p}}\bigg{)} \tag{7}\]  
\[\frac{dx_{p}}{dt}=v_{p} \tag{8}\]  
where the subscripts \(p\) and \(a\) denotes the particle and air, \(v\) is the velocity vector, \(g\) is the acceleration of gravity, \(x_{p}\) is the particle position, and the term \(\beta\big{(}v_{a}-v_{p}\big{)}\) is the drag of per unit particle mass, with \(\beta\) being calculated from:  
\[\beta=\frac{3\rho C_{D,p}|v_{a}-v_{p}|}{4\rho_{p}D_{p}} \tag{9}\]  
where \(D_{p}\) is the particle diameter and \(C_{D,p}\) is the particle drag coefficient, estimated by using the spherical drag law of [PERSON] [39].