# Aerodynamic Sensitivity Analysis for a Wind Turbine Airfoil in an Air-Particle Two-Phase Flow
## 2 Numerical Method for Air-Particle Two-Phase Flow Simulation  
The continuum air phase is simulated by solving the following incompressible N-S equations written in coordinate invariant formulation  
\[\begin{array}{c}\frac{\partial u_{i}}{\partial x_{i}}=0\\ \frac{\partial u_{i}}{\partial t}+u_{i}\frac{\partial u_{i}}{\partial x_{j}}=- \frac{1}{\rho}\frac{\partial p}{\partial x_{i}}+\
u\
abla^{2}u_{i}+\frac{1}{ \rho}Fp_{,i}\end{array} \tag{1}\]  
where \(x_{i}\) stands for a coordinate direction, \(u_{i}\) denotes a component of fluid velocity vector, \(\rho\), \(p\) and \(\
u\) are the fluid density, pressure and kinematic viscosity coefficient, respectively, \(\
abla^{2}\) represents the Laplace operator, and \(Fp_{,i}\) refers to the momentum source/sink \"body force\" term accounting for the particle effect on the fluid phase.  
The \(k-\omega\) Shear-Stress Transport (SST) two-equation model [19] is employed in the turbulence modeling, which combines the positive features of \(k-\omega\) model and a high Reynolds number \(k-\varepsilon\) model. The \(k-\omega\) SST model has been proven to work well for wind turbine applications [20]. The transport equations for the turbulent kinetic energy and the specific dissipation of turbulence read in differential form  
\[\begin{array}{c}\rho\frac{\partial k}{\partial t}+\rho\frac{\partial}{\partial x _{j}}\big{(}u_{j}k\big{)}=\frac{\partial}{\partial x_{j}}\Big{(}\Gamma_{k} \frac{\partial k}{\partial x_{j}}\Big{)}+G_{k}-Y_{k}\\ \rho\frac{\partial\omega}{\partial t}+\rho\frac{\partial}{\partial x_{j}}\Big{(} u_{j}\omega\Big{)}=\frac{\partial}{\partial x_{j}}\Big{(}\Gamma_{\omega} \frac{\partial\omega}{\partial x_{j}}\Big{)}+G_{\omega}-Y_{\omega}+D_{\omega} \end{array} \tag{2}\]  
where \(k\) is the turbulent kinetic energy, \(\omega\) is the rate of dissipation per unit turbulent kinetic energy. The terms on the right-hand side of Equation (2) represent conservative diffusion, eddy-viscosity production and dissipation, respectively, and the last term in the \(\omega\)-equation describes the cross diffusion. More details about the \(k-\omega\) SST model can be found in Ref. [19].  
For the continuum phase, the SIMPLEC (Semi-Implicit Method for Pressure-Linked Equations-Consistent) algorithm [21] and the second-order upwind scheme are used for the pressure-velocity coupling and the convection terms, respectively. An implicit method is applied to the time discretization.  
For the dispersed phase, 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 from the following DPM [13] formulated in a Lagrangian reference frame  
\[\begin{array}{c}\frac{dp_{P}}{dt}=\beta(\mathbf{v}-\mathbf{v}_{P})+g\Big{(}\frac{ \rho p-\rho}{\rho p}\Big{)}\\ \frac{dx_{P}}{dt}=\mathbf{v}_{P}\end{array} \tag{3}\]  
where the subscript \(P\) denotes the particle phase, \(\mathbf{v}\) is the velocity vector, \(\mathbf{g}\) is the acceleration of gravity, \(\mathbf{x}_{P}\) is the particle position, and the term \(\beta(\mathbf{v}-\mathbf{v}_{P})\) is the drag of per unit particle mass with \(\beta\) being calculated from  
\[\beta=\frac{3\rho C_{D,P}|\mathbf{v}-\mathbf{v}_{P}|}{4\rho P_{D}P} \tag{4}\]  
where \(D_{P}\) is the particle diameter, \(C_{D,P}\) is the particle drag coefficient estimated by using the spherical drag law of [PERSON] [22].  
Particle drag acts as the momentum coupling between the fluid and particulate phases. It is explicitly considered in the particle equation of motion. As for the continuous phase, the additional momentum source/sink term \(F_{P,j}\) in Equation (1) can be then evaluated from  
\[\mathbf{F}_{P}=-\frac{1}{\Omega_{cell}}\sum\beta(\mathbf{v}-\mathbf{v}_{P})m_{P} \tag{5}\]  
where \(\Omega_{cell}\) is the volume of the grid cell, \(m_{P}\) is the particle mass, and the sum is over all particles that traverse the cell.