## 2 Method
### Description of Flow Solvers
#### 2.3.2 CFD Fluent Solver  
The validation of solvers and mesh settings is needed before carrying out the optimization studies. The validation studies were performed at a Reynolds number of \(2.4\times 10^{6}\) and a Mach number of 0.10. The range of angles of attack was considered from \(0^{\circ}\) to \(16^{\circ}\) to include the stall angle. The aerodynamic coefficients and the pressure distribution were obtained and further used for solution validation. A C-shaped computational domain with structured grids, as seen in Figure 4, was used in domain discretization after the performance of the mesh convergence test. The turbulence model used in this study was the Transition \((\gamma-Re_{\theta})\) SST. This model uses a combination of SST K-\(\omega\) coupled with intermittency \(\gamma\) and transition onset Reynolds number. \(Re_{\theta}\) is the critical Reynolds number where the intermittency starts [78]. Four transport equations of the Transition \((\gamma-Re_{\theta})\) SST model are given below:  
\[\frac{\partial(\rho k)}{\partial t}+\frac{\partial\big{(}\rho U_{j}k\big{)}}{ \partial x_{j}}=P_{k}-D_{k}+\frac{\partial}{\partial x_{j}}\Bigg{[}(\mu+\sigma _{k}\mu_{t})\frac{\partial k}{\partial x_{j}}\Bigg{]} \tag{10}\]  
\[\frac{\partial(\rho\omega)}{\partial t}+\frac{\partial\big{(}\rho U_{j}\omega \big{)}}{\partial x_{j}}=P_{\omega}-D_{\omega}+\frac{\partial}{\partial x_{j }}\Bigg{[}(\mu+\sigma_{\omega}\mu_{t})\frac{\partial\omega}{\partial x_{j}} \Bigg{]}+2(1-F_{1})\frac{\rho\sigma_{\omega 2}}{\omega}\frac{\partial k}{ \partial x_{j}}\frac{\partial\omega}{\partial x_{j}} \tag{11}\]  
\[\frac{\partial(\rho\gamma)}{\partial t}+\frac{\partial\big{(}\rho U_{j} \gamma\big{)}}{\partial x_{j}}=P_{\gamma}-E_{\gamma}+\frac{\partial}{ \partial x_{j}}\Bigg{[}\Big{(}\mu+\frac{\mu_{t}}{\sigma_{\gamma}}\Big{)} \frac{\partial\gamma}{\partial x_{j}}\Bigg{]} \tag{12}\]  
\[\frac{\partial\big{(}\rho\overline{Re_{\theta t}}\big{)}}{\partial t}+\frac{ \partial\big{(}\rho U_{j}\overline{Re_{\theta t}}\big{)}}{\partial x_{j}}=P_{ \theta t}+\frac{\partial}{\partial x_{j}}\Bigg{[}\sigma_{\theta t}(\mu+\mu_{t}) \frac{\partial\overline{Re_{\theta t}}}{\partial x_{j}}\Bigg{]} \tag{13}\]The spatial discretization used in the solution is of the second order for all the turbulence statistical parameters, including the momentum, the turbulent kinetic energy, the specific dissipation rate, the intermittency and the Reynolds momentum thickness. ANSYS Fluent was used in this study to carry out the numerical analysis. The meshing was done using ANSYS Mesher, as shown in Figure 4. The rectangular computational domain was designed around the airfoil geometries. The outlet length was selected to be 30c, which was adequate to allow for the full development of the wake flow, based on the size of the computational domain adopted in previous studies. The distance from the domain inlet to the airfoil was 10c to prevent the inlet boundary from unphysically impacting the induction field upstream of the airfoil. The high-quality grids were generated with dense grids in the boundary layer of the airfoil with a gradual decrease of grid cells away from the airfoil surface. The first layer's thickness was calculated based on inlet velocity and adopted accordingly to capture the laminar and transitional boundary layer regions. Furthermore, the grid sensitivity analysis was performed according to the number of cells in the computational domain. A freestream velocity of 34 m/s with a turbulence rate of 0.01 percent was imposed at the inlet. The surface of the airfoil was modelled as a zero-roughness no-slip wall, and its values were executed on the grid domain, which prescribed its motion relative to the rest of the computational domain. A symmetry boundary condition was applied to achieve a parallel flow at the top and bottom of the domain by assuming zero normal velocity and zero normal gradients of the flow quantities. A zero static gauge pressure was applied at the exit of the domain.