Numerical Simulation of the Flow around NACA0018 Airfoil at High Incidences by Using RANS and DES Methods

Abstract

In this work, the flow around the NACA0018 airfoil with a wide range of attack angles was investigated based on the open-source computational fluid dynamics (CFD) platform OpenFOAM. Two numerical methods, Reynolds-averaged Navier–Stokes (RANS) and the detached eddy simulation (DES), were employed. Under the premise of a grid convergence analysis, the computed lift and drag coefficients were validated by the available experimental data. The pressure distribution, the complex flow mechanisms of the airfoil under the attached flow regime, the mild separation flow regime, and the post-stall flow regime, combined with the shedding vortex structures, streamlines, and vorticity distributions, are discussed. From the numerical results, it can be seen that the DES computation presents a better accuracy in the prediction of the lift and drag coefficients, with a deviation less than 10% at the largest angle of attack. Meanwhile, it also presents remarkable improvements in capturing the local flow field details, such as the unsteady separated flow and the shedding vortex structures.

1. Introduction

The issue of flow around airfoils at high incidences plays an important role in the research of ship and offshore structures, such as evaluating ship maneuverability and propulsion performance, optimizing wind turbine designs, etc. However, the flow field characteristics of airfoils at high incidences are rather complicated, due to the transition and pressure-induced separation. The flow features are highly related to the Reynolds number, the angle of attack, as well as the shape of the airfoil. Hence, it is important to explore the flow features around airfoils at high incidences.

Several research studies have been performed by using the experimental approach, aimed at exploring the flow mechanisms governing the complex flow around airfoils at high incidences. Jacobs and Sherman [1], Goett and Bullivant [2], and Sheldahl and Klimas [3] investigated the effect of the Reynolds number on the characteristics of an airfoil section. Nakano et al. [4] explored the flow separation and reattachment around an airfoil and the velocity field across the boundary layers over the airfoil surface. Timmer [5] measured the aerodynamic characteristics of an airfoil to serve as input to estimate the performance of vertical axis wind turbines. Boutilier et al. [6] studied the shear layer development over an airfoil using a combination of flow visualization, velocity field mapping, and surface pressure fluctuation measurements. Gim and Lee [7] studied the tip vortex structures around an airfoil with and without an endplate. Greenblatt et al. [8] examined pitching airfoils under relatively free stream oscillatory conditions (in a water tunnel and a wind tunnel, separately).

Concerning the high expenses and complex facilities of the experimental approach, the numerical approach, mainly referring to the computational fluid dynamics (CFD) method, has been widely adopted in this issue. Several research studies have used the Reynolds-averaged Navier–Stokes (RANS) method. Hassan et al. [9] investigated the accuracy of different turbulence models in predicting the aerodynamic characteristics of an airfoil. Raj [10] analyzed the flow over an airfoil at different Mach numbers and attack angles. Suvanjumrat [11] compared different turbulence models in the solution of the flow field around an airfoil. Yılmaz et al. [12] analyzed different airfoils at various attack angles and examined the effects of the symmetrical and asymmetrical profiles. Wang et al. [13] conducted a numerical study for a vertical axis wind turbine with different thicknesses or cambers of airfoils. Aqilah et al. [14] discussed the significant effect of mesh quality on the solution of the flow around an airfoil. Zidane et al. [15] explored the structure and behavior of the transitional separation bubbles on the airfoil at different Reynolds numbers. Li et al. [16] analyzed the effect of camber, the rotating axis position of the flap, the angle of attack, and flap thickness on the flow characteristics of a two-element wing sail. Overall, it was shown that the RANS method is suitable for the solution of attached flow, as well as a slightly separated flow, while it is hard to capture the large-scale flow separation for the airfoil at high incidences.

In contrast, the large eddy simulation (LES) method was proven to be capable of capturing large-scale flow separation. Breuer and Jovičić [17] investigated the separation flow past an airfoil at high attack angles. Mary and Sagaut [18] simulated the turbulent flow past an airfoil near the stall. Yuan et al. [19] investigated the flow past an airfoil, to clarify the viscous features of laminar separation and the transition flow followed by the complicated behavior of the flow reattachment. Li et al. [20] simulated the flow around an airfoil with high attack angles and discussed the vortex diffusion after flow separation. Breuer [21] investigated the effect of inflow turbulence on the flow around an airfoil with a laminar separation bubble. However, due to the high-computational resource consumption in the resolution of the near-wall flow, the LES method is too expensive for practical applications.

The hybrid RANS/LES method has been greatly developed in recent years to resolve the complex flow in practical applications (with affordable computational expenses). The method applies the RANS method for the attached flow in the boundary layer, and adopts the LES method for the separated flow regime, balancing the efficiency of the RANS method and the accuracy of the DES method. The widely used hybrid RANS/LES method and the detached eddy simulation (DES) method were proposed by Spalart et al. [22]. However, the original DES method revealed a defect of grid-spacing dependence; hence, it may behave incorrectly at the boundary layer and result in modeled stress depletion (MSD). To address the problem, two improved DES variants, delayed DES (DDES) [23] and improved DDES (IDDES) [24], were developed to ensure the transition of the RANS model to the LES mode was independent of grid spacing. Several research studies have been performed using the DES-type method. Schmidt and Thiele [25] investigated the influence of transient flow patterns on the quality of flow prediction. Li et al. [26] calculated pre-stall and post-stall aerodynamic characteristics of airfoils with different flow separation patterns at the stall regime. Probst et al. [27] and Im and Zha [28] discussed the effectiveness of different numerical schemes in the simulation of an airfoil at the stall. Grossi et al. [29] investigated the transonic flow over an airfoil within the buffet regime. Liang and Xue [30] investigated the wing-tip vortex from a rectangular wing with a square tip. Gan et al. [31] simulated the flow around a transonic wing flutter using a fully-coupled fluid/structure interaction (FSI) with high-order shock-capturing schemes. Xu et al. [32] investigated the flow over airfoils at a wide range of attacks. Yalçın et al. [33] explored the effects of two different length scale definitions in the simulation of the flow around an airfoil. Wang et al. [34] discussed the effects of the time step, spanwise lengths, and grid resolution on the prediction of flow around an airfoil beyond stall. Patel and Zha [35] simulated the post-stall flow around an airfoil to investigate the physics of flow separation. From the above research studies, it can be seen that the DES method is capable of predicting the unsteady flow with transition and separation.

Although large amounts of numerical research studies have been performed, the details and the corresponding discussions on the flow characteristics of airfoils under different regimes are relatively limited. Moreover, the effectiveness of the existing numerical approaches needs to be further explored, such as the difference between the RANS and DES methods in a complex flow solution with remarkable separation. In this paper, the flow around airfoils at high incidences was investigated based on the open-source CFD platform OpenFOAM. Using a symmetrical airfoil (NACA0018) as the study object, a series of numerical simulations were performed for a wide range of attack angles by using the RANS and DES methods. The computed lift and drag coefficients were compared with the available experimental data, to qualify the efficiency of the two numerical methods. Further, the flow features of the vortex structure, streamlines, and vorticity distribution, as well as the pressure distribution on the airfoil surface, are presented, to analyze the flow mechanism of the airfoil under the attached flow regime, the mild separation flow regime, and the post-stall flow regime. Based on the comparisons of the obtained flow field details, the capabilities of the RANS and DES solvers in OpenFOAM—in the prediction of the flow around airfoils at high incidences—were further assessed.

2. Numerical Method

For the adopted RANS method, the Spalart–Allmaras one-equation model was adopted for turbulence modeling. The adopted DES formulation was also derived from the Spalart–Allmaras one-equation model, and an improved variation of the original DES method, DDES method, was used in the study. In the Spalart–Allmaras one-equation model, the eddy viscosity $v_t$ is defined with a modified turbulence viscosity $\tilde{v}$, $v_t$ = $\tilde{v}{f}_{v1}$, and the transport equation is given by:

(1)$$\begin{gathered} \frac{D}{Dt}\left(\rho \tilde{v}\right)=\nabla ·\left(\rho D_{\tilde{v}}\tilde{v}\right)+\frac{C_{b2}}{{\sigma }_{v_t}}\rho {\left|\nabla \tilde{v}\right|}^2+ \\ {C}_{b1}\rho \tilde{S}\tilde{v}\left(1-f_{t2}\right)-\left(C_{w1}{f}_w-\frac{C_{b1}}{{\kappa }^2}{f}_{t2}\right)\rho \frac{{\tilde{v}}^2}{{\tilde{d}}^2}+S_{\tilde{v}} \end{gathered}$$

Here,

(2)$$\begin{gathered} f_{v1}=\frac{{\chi }^3}{{\chi }^3+C_{v1}^3}, \chi =\frac{\tilde{v}}{v} , C_{w1}=\frac{C_{b1}}{{\kappa }^2}+\frac{1+C_{b2}}{{\sigma }_{vt}} \\ {f}_{t2}=C_{t3}\exp \left(-C_{t4}{\chi }^2\right), f_w=g\left[\frac{1+C_{w3}^6}{g^6+C_{w3}^6}\right], g=r+C_{w2}\left(r^6-r\right) \\ r=\frac{\tilde{v}}{\tilde{S}{\kappa }^2{d}^2}, f_{v2}=1-\frac{\chi }{1+\chi f_{v1}}, \tilde{S}=\left|\omega \right|+\frac{\tilde{v}}{{\kappa }^2{d}^2}{f}_{v2} \end{gathered}$$

where $\rho$ is the fluid density, $d$ is the distance to the closest solid surface.

$\tilde{d}$ in the DDES model, different from the original DES model, is defined by:

(3)$\tilde{d}=\max \left[L_{\mathrm{RANS}}-f_d,\max \left(L_{\mathrm{RANS}}-L_{\mathrm{LES}}, 0\right)\right]$

The RANS length scale is defined by $L_{\mathrm{RANS}}$ = $d$, while the LES length scale is defined by $L_{\mathrm{LES}}$ = $\mathsf{\Psi }{C}_{DES}\Delta$, where $\Delta$ is the local grid scale, defined as $\Delta$ = $\sqrt[3]{\Delta x\times \Delta y\times \Delta z}$. $\mathsf{\Psi }$ is the low Reynolds number correction function:

(4)${\mathsf{\Psi }}^2=\min \left[{10}^2,\frac{1-\frac{1-C_{b1}}{C_{w1}{\kappa }^2{f}_w^{*}}\left[f_{t2}+\left(1-f_{t2}\right)f_{v2}\right]}{f_{v1}\max \left({10}^{-10},1-f_{t2}\right)}\right]$

The delay function is given by:

(5)$f_d=1-\tan \mathrm{h}\left[{\left(8r_d\right)}^3\right]$

Here, when $f_d$ = 0, the RANS model is recovered, and when $f_d$ = 1, the DES model is recovered. The $r_d$ parameter is given by:

(6)$r_d=\min \left(\frac{v+v_t}{\left|\nabla u\right|{\kappa }^2{d}^2},10\right)$

The model constants are given in Table 1; for more details on the SA-DDES model, refer to [23,36].

3. Case Description

3.1. Study Object and Working Condition

The symmetry airfoil NACA0018 was selected as the study object, and the profile is given in Figure 1. The chord length of $C$ = 0.154 m was adopted, which is in accordance with the published wind tunnel test [3]. Differently, the flowing medium was set to water with $\rho$ = 997.561 kg/m³ and $\mu$ = 8.89 × 10⁻⁴ Pa⋅s. The details of the studied working condition are summarized in Table 2.

For the sake of analysis, the non-dimensional lift and drag coefficients, $C_L$ and $C_D$, were adopted in the following, computed by:

(7)$C_L=L/\left(0.5\rho A{U}_0^2\right)$

(8)$C_D=D/\left(0.5\rho A{U}_0^2\right)$

where $L$ and $D$ denote the lift and drag acting on the airfoil. $A$ is the area of the airfoil, $A=C$×$S$. $U_0$ is the velocity of free inflow.

3.2. Computational Domain and Boundary Condition

A cuboid-shaped computational domain was adopted in the present numerical simulation, as shown in Figure 2. With the coordinate system defined at the center of the airfoil profile, the dimension of the computational domain is defined as: −4.9$C$ < $X$ < 9.7$C$, −4.9$C$ < $Y$ < 4.9$C$, and −0.5$S$ < $Z$ < 0.5$S$. As for the boundary condition, a uniform inflow condition with zero gradients for the pressure field was applied for the inlet side. A uniform outflow condition with the pressure field setting to the atmospheric pressure was applied for the outlet side. The slip wall boundary condition was set on the top and the bottom sides. The symmetry plane boundary condition was set on the front and the back sides. Moreover, the no-slip wall with zero velocity was set on the surface of the airfoil.

The discretization of the computational domain was performed by using the pre-processing tools in OpenFOAM. At first, the computational domain was discretized with a uniform hexahedral grid by using the blockMesh tool. Then, the grid refining, splitting, and snapping were performed by using the snappyHexMesh tool. A series of mesh refinement blocks were added to refine the grid spacing around the airfoil, especially the downstream of the airfoil. On the surface of the airfoil, five layers of prismatic cells were adopted to achieve a better resolution of the near-wall flow. The alignment of those cells was set in accordance with the targeted non-dimensional distance from the wall Y⁺ ranging from 30 to 60. See Figure 3 for the grid generation of the computational domain and the region around the airfoil.

3.3. Other Setups in the Numerical Solution

The present numerical study was performed by using the psioFoam solver in OpenFOAM. It is a pressure-based, coupled, unsteady solver, and solves the incompressible flow by using the pressure implicit with splitting of operator (PISO) algorithm. The discretization schemes adopted in the present study are given in Table 3, where $\varphi$ is the volumetric flux. For the temporal discretization, the first-order Euler scheme was used in the RANS computation, and the second backward scheme was used in the DES computation. As for the spatial discretization, the second-order Gauss integration scheme was used for divergence and gradient schemes in the RANS and DES computations, and the surface normal gradient scheme was adopted for the Laplacian scheme, which is a corrected unbounded second-order conservative scheme. Moreover, the first-order linear scheme was adopted for the interpolation scheme.

4. Convergence Study

A grid convergence study was performed with the methodology presented by Stern et al. [37] to clarify the effect of grid spacing in the numerical solution. Typically, the angle of attack $\alpha$ = 40 deg was selected as the study case. Three grids were generated based on the refinement factor of $\sqrt{2}$. For the RANS and DES computations, the time step was set to 1.0 × 10⁻³ s. To achieve the iterative convergence, the total simulated physical time was set to 5 s for the RANS computations and 10 s for the DES computations. A workstation with the CPU of Intel@ Xeon(R) Gold 5218 CPU@ 2.3 GHz (Intel Corp., Santa Clara, CA, USA) with 64 GB RAM was adopted in the present numerical study, and 16 processors were used for each computation. The computational details for the grid convergence study are given in Table 4, where the computational time refers to the wall-clock time.

The results of the grid convergence study are given in Table 5 and Table 6. According to the spatial discretization scheme, the theoretical order of accuracy $P_{Gest}$ = 2 was adopted. For the RANS computation, a relatively lower grid uncertainty was observed, with the values of $U_G\%S_{G1}$ ranging from 1.0% to 1.5% for $C_L$ and $C_D$. As for the DES computation, a lower grid uncertainty was presented for $C_L$ with the value of $U_G\%S_{G1}$ approaching 0.5%. However, a relatively higher grid uncertainty was presented for $C_D$, with the value of $U_G\%S_{G1}$ reaching 3.5%. The higher grid uncertainty in the DES computation was highly related to the larger values of ${\epsilon }_{32}\%S_{G1}$, resulting from the irrational solution with the coarse grid $S_{G3}$. It suggests that the coarse grid would be unable to solve the important flow features involved in the DES computation, as proposed by Muscari et al. [38]. A further refined grid would be helpful to improve the uncertainty, while it is rather time-consuming and unpractical. For another, the values of ${\epsilon }_{21}\%S_{G1}$ are at reasonably small magnitudes, demonstrating the minor difference between the solutions with the grid of $S_{G1}$ and $S_{G2}$. It can be concluded that the solution with the grid above the medium level would be less affected by the grid spacing. Therefore, the grid with the medium level was adopted in the following numerical studies, to balance the accuracy and efficiency of the computation.

5. Numerical Results

5.1. Lift and Drag Coefficients

Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11 plot the comparisons of the time histories of $C_L$ and $C_D$ obtained from the RANS and DES computations, and the last 2 s of the computations were selected. From the figures, it can be seen that $C_L$ and $C_D$, obtained from the RANS and DES computations, achieved stable values at small attack angles ($\alpha$ < 15 deg). At the angle of attack $\alpha$ = 15 deg, obvious fluctuations were presented for $C_L$ and $C_D$ obtained from the DES computation, which were not observed in the RANS computation. At the angle of attack $\alpha$ = 20 deg, both computations presented remarkably irregular fluctuations, and the frequency of fluctuations was relatively reduced for the DES computation. As the angle of attack increased further, $C_L$ and $C_D$, obtained from the RANS computation, tended to present regular fluctuations, and the amplitudes of the fluctuations were enlarged with the increase of the attack angles. Differently, the tendencies of $C_L$ and $C_D,$ obtained from the DES computations, were still chaotic, and the amplitudes were relatively less affected by the attack angles.

Figure 12 plots the comparison of the time-averaged $C_L$ and $C_D$ obtained from the RANS and DES computations. The numerical results are also compared with the available experimental data [3,5]. The time-averaged values of $C_L$ and $C_D$ take the averaged values plotted in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11. From the experimental data, it can be seen that $C_L$ presents a monotonically increasing tendency at $\alpha$ < 15 deg, and the tendency is approximately linear. After $\alpha$ > 15 deg, there is a remarkable decline/fluctuant in the tendency. Differently, $C_D$ presents a small value at the small attack angles ($\alpha$ < 15 deg), and it is enlarged remarkably after $\alpha$ > 15 deg. Hence, the angle of attack $\alpha$ = 15 deg is near the critical point of the airfoil stall, and the condition of $\alpha$ > 15 deg is known as the post-stall regime. Moreover, it is shown that a remarkable difference is observed between the experimental data from Sheldahl and Klimas [3] and Timmer [5]. In particular, a stall hysteresis loop was observed in the experimental data by Timmer [5]. The differences may derive from the differences in test conditions, while limited information is published. It is shown that the RANS and DES computations present rather different tendencies of $C_L$ and $C_D$ with the variations of the attack angles, especially for the post-stall regime.

The deviations of $C_L$ and $C_D$ between the numerical results and the experimental data are given in Table 7 and Table 8. By comparing the numerical results and the experimental data, it can be seen that for the smaller attack angles ($\alpha$ ≤ 15 deg), the RANS computation presents better consistency with the experimental data and the deviation of $C_L$ is less than 5% for most attack angles. In contrast, the DES computation presents a higher deviation, with the error reaching about 20%, especially for the angle of attack $\alpha$ = 15 deg. Moreover, the deviation of $C_D$ is relatively higher due to the small magnitude before the stall, reaching about 150% for the RANS computation and about 200% for the DES computation. As the angle of attack increases further, the deviations of $C_L$ and $C_D$ between the RANS computation and the experimental data are enlarged, reaching about 30% at the largest angle of attack $\alpha$ = 40 deg. Differently, the deviations between the DES computation and the experimental data are reduced, with the values reaching about 1% for $C_L$ and about 10% for $C_D$ at the angle of attack $\alpha$ = 40 deg. Based on the comparison between the numerical results and the experimental data, it can be concluded that the RANS computation has higher accuracy in the prediction of the global loads before the airfoil stall, as well as the critical point of the airfoil stall, while the DES computation presents a higher accuracy in the prediction of the global loads after the airfoil stall.

5.2. Shedding Vortex Structures

Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18, Figure 19 and Figure 20 present the comparisons of the transient vortex structures shedding from the airfoil obtained from the RANS and DES computations. The vortex structures are visualized with the isosurface of the second invariant of the rate of strain tensor $Q$, and the isosurface is colored by the non-dimensional axial velocity $U_x/U_0$. From Figure 13 and Figure 14, it can be seen that there are almost no vortex structures generated from the airfoil surface under the small attack angles ($\alpha$ < 15 deg), and the near-wall flow is reduced along the airfoil surface. Moreover, a minor difference is observed for the RANS and DES computations. As shown in Figure 15, a slight disturbance of the near-wall flow appears in the RANS computation at the angle of attack $\alpha$ = 15 deg, but a rather weak shedding vortex is presented. Differently, a series of remarkable vortices were generated from the latter part of the airfoil in the DES computation, and the vortex structures developed well along the downstream.

From Figure 16, Figure 17, Figure 18, Figure 19 and Figure 20, one can see that large-scale shedding vortices are presented in the RANS and DES computations, as the angle of attack increases further ($\alpha$ > 15 deg). The vortex structures present remarkable differences between the two numerical methods. In the RANS computation, two large-scale shedding vortices were generated near the leading and trailing edges of the airfoil, respectively. The vortex structures are relatively simple, and obvious spanwise uniformity is observed. By contrast, more vortices in different scales are observed in the DES computation, and the locations of the shedding vortices cover the whole upper surface of the airfoil. The vortex structures are rather complicated and present remarkable spanwise non-uniformity. In particular, the breakup of the large-scale vortices is well captured by the DES computation, which is predicted unrealistically in the RANS computation.

In combination with the computed lift and drag coefficients, it can be concluded that the variations of the global loads on the airfoil are highly related to the generation and development of the shedding vortices. For the small angle of attack ($\alpha$ < 15 deg), there is no obvious vortex structure around the airfoil. Correspondingly, the drag of the airfoil has a rather small magnitude. At the angle of attack $\alpha$ > 15 deg, remarkable vortex structures appear and result in a significant increase of the drag. Meanwhile, due to the generation and development of the shedding vortices, the lift and drag coefficients present obvious fluctuation tendencies. The fluctuations are regular for the RANS computation and are chaotic for the DES computation, corresponding to the developments of the shedding vortex structures, respectively.

5.3. Streamlines and Vorticity Distribution

Figure 21, Figure 22, Figure 23, Figure 24, Figure 25, Figure 26, Figure 27 and Figure 28 present the comparisons of the transient streamlines and vorticity contours around the airfoil profiles obtained from the RANS and DES computations. The mid-span section of $Z$ = 0 was taken into account. Based on the streamlines and vorticity contours, the attached flow regime, mild separation flow regime, and post-stall flow regime were identified clearly. As shown in Figure 21 and Figure 22, the attached flow was the dominant feature for the flow around the airfoil at the small angle of attack ($\alpha$ < 15 deg). It is shown that the local flow is smooth over the surface of the airfoil and no flow separation occurs. Moreover, the flow feature is consistent for the RANS and DES computations.

As the angle of attack increases to $\alpha$ = 15 deg, the attached flow regime transforms into the mild separation flow regime. From Figure 23, one can see that mild flow separation occurs and an obvious separation bubble is formed near the trailing edge. Meanwhile, an apparent difference is observed between the flow features captured by the RANS and DES computations. For the RANS computation, the flow separation point is near the mid-chord of the airfoil. The separation bubble is formed by a clockwise vortex shedding from the separation point, and the local flow around the separation bubble is relatively smooth. As for the DES computation, the separation point is located slightly ahead of that in the RANS computation. The separation bubble has a larger scale, which is formed by a couple of vortices, including a clockwise vortex shedding from the separation point and an anticlockwise vortex generated near the trailing edge.

At the larger angle of attack ($\alpha$ > 15 deg), the mild separation flow regime transformed into the post-stall regime. It is shown that the flow separation on the upper surface of the airfoil was intensified, and the distribution of separated flow changed remarkably with the increase of the angle of attack. Remarkable differences appeared between the flow fields from the RANS and DES computations. For the RANS computation, the separation point moved forward to the leading edge at the attack angles $\alpha$ = 20 deg and 25 deg. Consequently, a separation bubble covering the whole upper surface was formed, including a clockwise vortex shedding from the separation point and an anticlockwise vortex generated near the trailing edge. As the angle of attack increased further, another derived shedding vortex appeared near the leading edge. Overall, the local flow fields around the airfoil were relatively regular in the RANS computation, which was dominated by vortices shedding from the separation point and the trailing edge. In contrast, the local flow presented by the DES computation was rather complicated. The flow separation region covered the whole upper surface of the airfoil, and the separation bubble was formed by a series of chaotic detached vortices. Different from the regular diffusion of vortices presented by the RANS computation, the breakup of the large-scale shedding vortices was captured well by the DES computation. From the streamlines and vorticity distributions, it can be seen that, due to the generation and development of the shedding vortices, as well as the interaction between the vortices in different directions, strong unsteady characteristics were observed for the local flow around the airfoil at the larger angle of attack.

5.4. Pressure Distribution

Figure 29 presents the comparisons of the transient pressure distributions over the airfoil surface obtained from the RANS and DES computations. Corresponding to the streamlines and vorticity distributions discussed above, the mid-span section of $Z$ = 0 was taken into account. The pressure on the surface of the airfoil is represented with the non-dimensional pressure coefficient $C_p=p/\left(0.5\rho U_0^2\right)$, where $p$ is the relative pressure. From Figure 29a,b one can see that the discrepancies of the pressure distributions between the RANS and DES computations were rather limited at the small attack angles ($\alpha$ < 15 deg), coinciding with the obtained streamlines and vorticity distribution.

As the angle of attack increased to $\alpha$ = 15 deg, an obvious difference in the pressure distribution between the RANS and DES computations was observed. According to the above analysis of the flow field details, a mild flow separation appeared under the condition, and the separation points in the RANS and DES computations were located at about $x/C$ = 0.5. From Figure 29c, it can be seen that the negative pressure on the upper surface was reduced after the separated point. Meanwhile, due to the small recirculation flow near the trailing edge (as shown in Figure 23), a small negative pressure region appeared on the upper and lower surfaces near the trailing edge, and it was more distinct for the DES computation. Overall, the DES computation presents a weaker negative pressure on a major part of the upper surface, resulting in a lower lift than that of the RANS computation.

For the attack angles $\alpha$ = 20 deg and 25 deg, the difference in the pressure distribution on the airfoil surface between the RANS and DES computations was more obvious, especially for the upper surface near the leading edge. Under the conditions, the separation point of the RANS computation was about $x/C$ = 0.15 for $\alpha$ = 20 deg and about $x/C$ = 0.05 for $\alpha$ = 25 deg, while the separation point of the DES computation was approaching the leading edge. From Figure 29d,e it can be seen that the negative pressure on the upper surface was reduced and presented a flat distribution for the RANS computation after the separated point. Differently, the negative pressure on the upper surface was relatively irregular, in accordance with the chaotic local flow field. Moreover, the DES computation presented a stronger negative pressure on the upper surface. As a result, the lift obtained from the DES computation was higher than that of the RANS computation. For the larger attack angles $\alpha$ = 30 deg, 35 deg, and 40 deg, the incorrectly separated flow and the diffused vortices in the RANS computation resulted in a stronger negative pressure on the upper surface of the airfoil, which consequently led to a much higher lift and drag coefficients. In contrast, the pressure distribution obtained from the DES computation was more reasonable, according to the comparison of the lift and drag coefficients with the experimental data.

6. Conclusions

In this paper, the flow around the airfoil NACA0018 at high incidences was investigated by using the open-source CFD platform OpenFOAM. RANS and DES methods were adopted to simulate the local flow fields around the airfoil, with the attack angles covering the attached flow regime, mild separation flow regime, and post-stall flow regime. Under the premise of a grid convergence analysis, the effectiveness of the numerical methods was discussed based on the comparison between the numerical results and the available experimental data. It is shown that the feasibility of the numerical methods highly depends on the flow regime. For the attached flow regime and mild separation flow regime, $C_L$ and $C_D$ obtained from the RANS computation presented an agreement with the experimental data, with the deviation less than 5% for most attack angles. As for the post-stall flow regime, the DES computation had high accuracy, with the deviation reaching about 1% for $C_L$ and about 10% for $C_D$ at the largest angle of attack $\alpha$ = 40 deg.

We discuss the obtained flow features of the shedding vortex structures, streamlines, and vorticity distribution, as well as the pressure distribution, the flow mechanism of the airfoil at high incidences, and the effectiveness of the adopted numerical methods. It is shown that vortex structures are relatively simple and present obvious spanwise uniformity in the RANS computation. In contrast, vortex structures with different scales are captured by the DES computation, and the vortex structures present remarkable spanwise non-uniformity. In particular, different from the regular diffusion of vortices presented by the RANS computation, the DES computation reproduces the breakup and development of large-scale shedding vortices. Based on the pressure distribution on the airfoil section, the feature of the flow separation and its effect on the global loads were explored. Based on the pressure distribution at high incidences, the unrealistic flow separation and shedding vortices in the RANS computation resulted in a too-strong suction pressure on the upper surface of the airfoil. Hence, $C_L$ and $C_D$ were remarkably overestimated. Differently, the DES computation presented a more reasonable pressure distribution on the airfoil surface. As a result, a good agreement was achieved between the DES computation and the experimental data. From the numerical results, it can be seen that the adopted DES method is better capable of predicting the flow around the airfoil at high incidences, especially for the post-stall regime, with remarkable flow separation and large-scale shedding vortices.

This study focused on the effectiveness of RANS and DES methods in the numerical solution of the flow around an airfoil at high incidences. Although the adopted DES method presented satisfactory accuracy overall, an obvious deviation appeared near the critical point of the airfoil stall. In the future, more systematic studies should be performed, exploring the suitable turbulence model and the DES variant, to further improve the prediction accuracy of the flow around the airfoil at high incidences.

Author Contributions

Conceptualization, H.G., G.L. and Z.Z.; methodology, H.G.; software, H.G. and G.L.; validation, H.G., G.L. and Z.Z.; formal analysis, H.G. and Z.Z.; investigation, H.G. and G.L.; data curation, H.G. and G.L.; writing—original draft preparation, H.G., G.L. and Z.Z.; writing—review and editing, H.G. and Z.Z.; supervision, Z.Z.; project administration, G.L. and Z.Z.; funding acquisition, H.G., G.L. and Z.Z. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Profile of the NACA0018 airfoil.

Figure 2. Illustration of the computational domain and boundary conditions.

Figure 3. Grid generation of (a) the computational domain and (b) the region around the airfoil.

Figure 4. Time histories of (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 5 deg.

Figure 5. Time histories of (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 10 deg.

Figure 6. Time histories of (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 15 deg.

Figure 7. Time histories of (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 20 deg.

Figure 8. Time histories of (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 25 deg.

Figure 9. Time histories of (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 30 deg.

Figure 10. Time histories of (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 35 deg.

Figure 11. Time histories of (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 40 deg.

Figure 12. Comparison of (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.] between the numerical results and the experimental data [3,5].

Figure 12. Comparison of (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.] between the numerical results and the experimental data [3,5].

Figure 13. [Forumla omitted. See PDF.] = 10,000 of the airfoil with [Forumla omitted. See PDF.] = 5 deg: (a) RANS; (b) DES.

Figure 14. [Forumla omitted. See PDF.] = 10,000 of the airfoil with [Forumla omitted. See PDF.] = 10 deg: (a) RANS; (b) DES.

Figure 15. [Forumla omitted. See PDF.] = 10,000 of the airfoil with [Forumla omitted. See PDF.] = 15 deg: (a) RANS; (b) DES.

Figure 16. [Forumla omitted. See PDF.] = 10,000 of the airfoil with [Forumla omitted. See PDF.] = 20 deg: (a) RANS; (b) DES.

Figure 17. [Forumla omitted. See PDF.] = 10,000 of the airfoil with [Forumla omitted. See PDF.] = 25 deg: (a) RANS; (b) DES.

Figure 18. [Forumla omitted. See PDF.] = 10,000 of the airfoil with [Forumla omitted. See PDF.] = 30 deg: (a) RANS; (b) DES.

Figure 19. [Forumla omitted. See PDF.] = 10,000 of the airfoil with [Forumla omitted. See PDF.] = 35 deg: (a) RANS; (b) DES.

Figure 20. [Forumla omitted. See PDF.] = 10,000 of the airfoil with [Forumla omitted. See PDF.] = 40 deg: (a) RANS; (b) DES.

Figure 21. Streamlines on the mid-span section with [Forumla omitted. See PDF.] = 5 deg: (a) RANS; (b) DES.

Figure 22. Streamlines on the mid-span section with [Forumla omitted. See PDF.] = 10 deg: (a) RANS; (b) DES.

Figure 23. Streamlines on the mid-span section with [Forumla omitted. See PDF.] = 15 deg: (a) RANS; (b) DES.

Figure 24. Streamlines on the mid-span section with [Forumla omitted. See PDF.] = 20 deg: (a) RANS; (b) DES.

Figure 25. Streamlines on the mid-span section with [Forumla omitted. See PDF.] = 25 deg: (a) RANS; (b) DES.

Figure 26. Streamlines on the mid-span section with [Forumla omitted. See PDF.] = 30 deg: (a) RANS; (b) DES.

Figure 27. Streamlines on the mid-span section with [Forumla omitted. See PDF.] = 35 deg: (a) RANS; (b) DES.

Figure 28. Streamlines on the mid-span section with [Forumla omitted. See PDF.] = 40 deg: (a) RANS; (b) DES.

Figure 29. Pressure distribution on the mid-span section of the airfoil at (a) [Forumla omitted. See PDF.] = 5 deg, (b) [Forumla omitted. See PDF.] = 10 deg, (c) [Forumla omitted. See PDF.] = 15 deg, (d) [Forumla omitted. See PDF.] = 20 deg, (e) [Forumla omitted. See PDF.] = 25 deg, (f) [Forumla omitted. See PDF.] = 30 deg, (g) [Forumla omitted. See PDF.] = 35 deg, and (h) [Forumla omitted. See PDF.] = 40 deg.

Table 1

Constants in the SA-DDES model.

${\mathit{\sigma }}_{\mathit{v}\mathit{t}}$	${\mathit{C}}_{\mathit{b}1}$	${\mathit{C}}_{\mathit{b}2}$	${\mathit{C}}_{\mathit{w}2}$	${\mathit{C}}_{\mathit{w}3}$	${\mathit{C}}_{\mathit{v}1}$	${\mathit{C}}_{\mathit{s}}$	${\mathit{C}}_{\mathit{D}\mathit{E}\mathit{S}}$	${\mathit{C}}_{\mathit{k}}$	${\mathit{C}}_{\mathit{t}3}$	${\mathit{C}}_{\mathit{t}4}$	${\mathit{f}}_{\mathit{w}}^{*}$
2/3	0.1355	0.622	0.3	2.0	7.1	0.3	0.65	0.07	1.2	0.5	0.424

Table 2

Study object and working conditions.

Terms	Values
Airfoil profile	NACA0018
Chord length $C$ (m)	0.154
Spanwise length $S$ (m)	0.1
Reynolds number $Re$	1,000,000
Angle of attack $\alpha$ (deg)	5, 10, 15, 20, 25, 30, 35, 40

Table 3

Discretization schemes adopted in the present study.

Terms		RANS	DES
Time		Euler	Backward
Divergence	$\left(\varphi ,U\right)$	Gauss linearUpwind grad(U)	Gauss LUST unlimitedGrad(U)
	$\left(\varphi ,\tilde{v}\right)$	Gauss limitedLinear 1	Gauss limitedLinear 1
Gradient	$U$	Gauss linear	cellLimited Gauss linear 1
	$\tilde{v}$	Gauss linear	cellLimited Gauss linear 1
Laplacian		Gauss linear limited corrected 0.5	Gauss linear limited corrected 0.33
Interpolation		linear	linear

Table 4

Computational details for the grid convergence study.

	Method	Fine Grid ${\mathit{S}}_{\mathit{G}1}$	Medium Grid ${\mathit{S}}_{\mathit{G}2}$	Coarse Grid ${\mathit{S}}_{\mathit{G}3}$
Grid number	-	3.98 × 106	2.07 × 106	1.06 × 106
Computationalexpense (h)	RANS	50–55	25–30	5–10
	DES	60–65	40–45	15–20

Table 5

Grid convergence analysis of the RANS computation.

	${\mathit{S}}_{\mathit{G}1}$	${\mathit{S}}_{\mathit{G}2}$	${\mathit{S}}_{\mathit{G}3}$	${\mathit{\epsilon }}_{21}\mathit{\%}{\mathit{S}}_{\mathit{G}1}$	${\mathit{\epsilon }}_{32}\mathit{\%}{\mathit{S}}_{\mathit{G}1}$	${\mathit{R}}_{\mathit{G}}$	${\mathit{P}}_{\mathit{G}}$	${\mathit{C}}_{\mathit{G}}$	${\mathit{U}}_{\mathit{G}}\mathit{\%}{\mathit{S}}_{\mathit{G}1}$
$C_L$	1.334	1.318	1.349	−1.210	2.303	−0.525	-	-	1.152
$C_D$	1.156	1.151	1.186	−0.373	2.981	−0.125	-	-	1.491

Table 6

Grid convergence analysis of the DES computation.

	${\mathit{S}}_{\mathit{G}1}$	${\mathit{S}}_{\mathit{G}2}$	${\mathit{S}}_{\mathit{G}3}$	${\mathit{\epsilon }}_{21}\mathit{\%}{\mathit{S}}_{\mathit{G}1}$	${\mathit{\epsilon }}_{32}\mathit{\%}{\mathit{S}}_{\mathit{G}1}$	${\mathit{R}}_{\mathit{G}}$	${\mathit{P}}_{\mathit{G}}$	${\mathit{C}}_{\mathit{G}}$	${\mathit{U}}_{\mathit{G}}\mathit{\%}{\mathit{S}}_{\mathit{G}1}$
$C_L$	1.015	1.017	1.116	0.222	9.709	0.023	10.898	42.684	0.439
$C_D$	0.802	0.819	0.886	2.095	8.468	0.247	4.030	3.042	3.501

Table 7

Deviation of $C_L$ between the numerical results and experimental data.

$\mathit{\alpha }$ (deg).	RANS	ERR%(Sheldahl)	ERR%(Timmer)	DES	ERR%(Sheldahl)	ERR%(Tmmer)
5	0.406	26.257	16.354	0.412	25.107	15.049
10	0.959	1.648	1.273	0.895	8.252	7.902
15	1.048	1.378	4.855	0.819	22.948	25.665
20	0.783	19.328	22.523	1.122	−15.557	−10.981
25	0.861	14.839	11.534	1.058	−4.729	−8.794
30	1.162	−35.867	-	1.073	−25.457	-
35	1.206	−23.025	-	1.024	−4.469	-
40	1.320	−27.559	-	1.020	1.469	-

Table 8

Deviation of $C_D$ between the numerical results and experimental data.

$\mathit{\alpha }$ (deg).	RANS	ERR%(Sheldahl)	ERR%(Timmer)	DES	ERR%(Sheldahl)	ERR%(Tmmer)
5	0.016	−73.263	−63.781	0.015	−60.580	−51.792
10	0.024	−56.362	−14.068	0.023	−46.694	−7.015
15	0.061	−142.182	4.196	0.079	−213.138	−23.873
20	0.249	11.695	−72.966	0.313	−11.055	−117.528
25	0.464	−14.461	−93.758	0.412	−1.829	−72.376
30	0.752	−31.855	-	0.536	5.973	-
35	0.914	−22.661	-	0.668	10.357	-
40	1.157	−25.798	-	0.817	11.179	-

References

1. Jacobs, E.N.; Sherman, A. Airfoil Section Characteristics as Affected by Variations of the Reynolds Number; Annual Report-National Advisory Committee for Aeronautics Iowa State University: Ames, IA, USA, 1937; 586.

2. Goett, H.J.; Bullivant, W.K. Tests of NACA 009, 0012 and 0018 Airfoils in the Full-Scale Tunnel; Annual Report-National Advisory Committee for Aeronautics Iowa State University: Ames, IA, USA, 1938; 647.

3. Sheldahl, R.E.; Klimas, P.C. Aerodynamic Characteristics of Seven Symmetrical Airfoil Sections through 180-Degree Angle of Attack for Use in Aerodynamic Analysis of Vertical Axis Wind Turbines; Sandia National Laboratories Energy Report Sandia National Labs.: Albuquerque, NM, USA, 1981.

4. Nakano, T.; Fujisawa, N.; Oguma, Y.; Takagia, Y.; Leeb, S. Experimental study on flow and noise characteristics of NACA0018 airfoil. J. Wind Eng. Ind. Aerod.; 2007; 95, pp. 511-531. [DOI: https://dx.doi.org/10.1016/j.jweia.2006.11.002]

5. Timmer, W.A. Two-dimensional low-Reynolds number wind tunnel results for airfoil NACA 0018. Wind Eng.; 2008; 32, pp. 525-537. [DOI: https://dx.doi.org/10.1260/030952408787548848]

6. Boutilier, M.S.H.; Yarusevych, S. Separated shear layer transition over an airfoil at a low Reynolds number. Phys. Fluids; 2012; 24, 084105. [DOI: https://dx.doi.org/10.1063/1.4744989]

7. Gim, O.S.; Lee, G.W. Flow characteristics and tip vortex formation around a NACA 0018 foil with anendplate. Ocean Eng.; 2013; 60, pp. 28-38. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2012.12.009]

8. Greenblatt, D.; Mueller-Vahl, H.; Strangfeld, C.; Medina, A.; Ol, M.; Granlund, K. High advance-ratio airfoil streamwise oscillations: Wind tunnel vs. water tunnel. Proceedings of the 54th AIAA Aerospace Sciences Meeting; San Diego, CA, USA, 4–8 January 2016.

9. Hassan, G.E.; Hassan, A.; Youssef, M.E. Numerical investigation of medium range re number aerodynamics characteristics for NACA0018 airfoil. CFD Lett.; 2014; 6, pp. 175-187.

10. Raj, J. CFD analysis of flow characteristics of NACA0012 airfoil using SU2. J. Mech. Aeronaut. Eng. Res.; 2017; 1, pp. 25-28.

11. Suvanjumrat, C. Comparison of turbulence models for flow past NACA0015 airfoil using OpenFOAM. Eng. J.; 2017; 21, pp. 207-221. [DOI: https://dx.doi.org/10.4186/ej.2017.21.3.207]

12. Yilmaz, M.; Koten, H.; Çetinkaya, E.; Coşar, Z. A comparative CFD analysis of NACA0012 and NACA4412 airfoils. J. Energy Syst.; 2018; 2, pp. 145-159. [DOI: https://dx.doi.org/10.30521/jes.454193]

13. Wang, Y.; Shen, S.; Li, G.; Huang, D.; Zheng, Z. Investigation on aerodynamic performance of vertical axis wind turbine with different series airfoil shapes. Renew. Energ.; 2018; 126, pp. 801-818. [DOI: https://dx.doi.org/10.1016/j.renene.2018.02.095]

14. Aqilah, F.; Islam, M.; Juretic, F.; Guerrero, J.; Wood, D.; Nasir, F. Study of mesh quality improvement for CFD analysis of an airfoil. IIUM Eng. J.; 2018; 19, pp. 203-212. [DOI: https://dx.doi.org/10.31436/iiumej.v19i2.905]

15. Zidane, I.F.; Swadener, G.; Saqr, K.M.; Ma, X.; Shehade, M.F. CFD investigation of transitional separation bubble characteristics on NACA 63415 airfoil at low Reynolds numbers. Proceedings of the 25th UKACM Conference on Computational Mechanics; Birmingham, UK, 12–13 April 2017.

16. Li, C.; Wang, H.; Sun, P. Numerical investigation of a two-element wingsail for ship auxiliary propulsion. J. Mar. Sci. Eng.; 2020; 8, 333. [DOI: https://dx.doi.org/10.3390/jmse8050333]

17. Breuer, M.; Jovičić, N. An LES investigation of the separated flow past an airfoil at high angle of attack. Direct and Large-Eddy Simulation IV; Springer: Dordrecht, The Netherlands, 2001; pp. 165-172.

18. Mary, I.; Sagaut, P. Large eddy simulation of flow around an airfoil near stall. AIAA J.; 2002; 40, pp. 1139-1145. [DOI: https://dx.doi.org/10.2514/2.1763]

19. Yuan, W.; Xu, H.; Khalid, M.; Radespiel, R. A parametric study of LES on laminar-turbulent transitional flows past an airfoil. Int. J. Comut. Fluid Dyn.; 2006; 20, pp. 45-54. [DOI: https://dx.doi.org/10.1080/10618560600578492]

20. Li, C.; Zhu, S.; Xu, Y.; Xiao, Y. 2.5 D large eddy simulation of vertical axis wind turbine in consideration of high angle of attack flow. Renew. Energ.; 2013; 51, pp. 317-330. [DOI: https://dx.doi.org/10.1016/j.renene.2012.09.011]

21. Breuer, M. Effect of inflow turbulence on an airfoil flow with laminar separation bubble: An LES study. Flow Turbul. Combust.; 2018; 101, pp. 433-456. [DOI: https://dx.doi.org/10.1007/s10494-017-9890-2]

22. Spalart, P.R. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach. Proceedings of the 1st AFOSR International Conference on DNS/LES; Ruston, LA, USA, 4–8 August 1997.

23. Spalart, P.R.; Deck, S.; Shur, M.; Squires, K.D. A new version of detached-eddy simulation, resistant to ambiguous grid densities. Theor. Comput. Fluid Dyn.; 2006; 20, pp. 181-195. [DOI: https://dx.doi.org/10.1007/s00162-006-0015-0]

24. Shur, M.L.; Spalart, P.R.; Strelets, M.K.; Travin, A.K. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities. Int. J. Heat. Fluid Flow; 2008; 29, pp. 1638-1649. [DOI: https://dx.doi.org/10.1016/j.ijheatfluidflow.2008.07.001]

25. Schmidt, S.; Thiele, F. Detached eddy simulation of flow around A-airfoil. Flow Turbul. Combust.; 2003; 71, pp. 261-278. [DOI: https://dx.doi.org/10.1023/B:APPL.0000014933.66058.22]

26. Li, D.; Men’shov, I.; Nakamura, Y. Detached-eddy simulation of three airfoils with different stall onset mechanisms. J. Aircr.; 2006; 43, pp. 1014-1021. [DOI: https://dx.doi.org/10.2514/1.17274]

27. Probst, A.; Wolf, C.; Radespiel, R.; Knopp, T.; Schwamborn, D. A comparison of detached-eddy simulation and Reynolds-stress modeling applied to the flow over a backward-facing step and an airfoil at stall. Proceedings of the 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition; Orlando, FL, USA, 4–7 January 2010.

28. Im, H.; Zha, G. Delayed detached eddy simulation of a stall flow over NACA0012 airfoil using high order schemes. Proceedings of the 49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition; Orlando, FL, USA, 4–7 January 2011.

29. Grossi, F.; Braza, M.; Hoarau, Y. Delayed detached-eddy simulation of the transonic flow around a supercritical airfoil in the buffet regime. Progress in Hybrid RANS-LES Modelling; Springer: Berlin/Heidelberg, Germany, 2012; pp. 369-378.

30. Liang, Z.; Xue, L. Detached-eddy simulation of wing-tip vortex in the near field of NACA 0015 airfoil. J. Hydrodynam. B; 2014; 26, pp. 199-206. [DOI: https://dx.doi.org/10.1016/S1001-6058(14)60022-6]

31. Gan, J.; Im, H.S.; Chen, X.; Zha, G.C.; Pasiliao, C.L. Delayed detached eddy simulation of wing flutter boundary using high order schemes. J. Fluids. Struct.; 2017; 71, pp. 199-216. [DOI: https://dx.doi.org/10.1016/j.jfluidstructs.2017.03.002]

32. Xu, H.Y.; Qiao, C.L.; Yang, H.Q.; Ye, Z.Y. Delayed detached eddy simulation of the wind turbine airfoil S809 for angles of attack up to 90 degrees. Energy; 2017; 118, pp. 1090-1109. [DOI: https://dx.doi.org/10.1016/j.energy.2016.10.131]

33. Yalçın, Ö.; Cengiz, K.; Özyörük, Y. High-order detached eddy simulation of unsteady flow around NREL S826 airfoil. J. Wind. Eng. Ind. Aerodyn.; 2018; 179, pp. 125-134. [DOI: https://dx.doi.org/10.1016/j.jweia.2018.05.017]

34. Wang, Y.; Liu, K.; Han, Z. Detached eddy simulation on the flow around NACA0021 airfoil beyond stall using OpenFOAM. Proceedings of the Asia-Pacific International Symposium on Aerospace Technology; Chengdu, China, 16–18 October 2018.

35. Patel, P.; Zha, G. Improved delayed detached eddy simulation of separated flow. Proceedings of the AIAA Aviation 2020 Forum; Virtual Event, 15–19 June 2020.

36. OpenFOAM. Available online: https://www.openfoam.com (accessed on 24 August 2021).

37. Stern, F.; Wilson, R.; Shao, J. Quantitative V&V of CFD simulations and certification of CFD codes. Int. J. Numer. Methods Fluids.; 2006; 50, pp. 1335-1355.

38. Muscari, R.; Di Mascio, A.; Verzicco, R. Modeling of vortex dynamics in the wake of a marine propeller. Comput. Fluids.; 2013; 73, pp. 65-79. [DOI: https://dx.doi.org/10.1016/j.compfluid.2012.12.003]