Development of a new aerofoil profile with a high lift-to-drag ratio for wind turbines using a low fidelity accurate optimization flow solver

Introduction

Owing to the increase in population and urbanization in developing countries, electric power requirements have raised significantly all over the world. The dependence of power generation on depleting fossil fuels and natural gas resources side by side with ever-growing environmental concerns demands exploration of renewable and environment-friendly power generation methods (WEC 2013). Renewable energy methods involving hydroelectric, solar power conversion, and wind power extraction are being actively carried out worldwide to help meet the ever-growing electric power requirements (Chakraborty and Razzak 2014).

The aerodynamic design of the wind turbine blades is one of the principal problems of the wind industry, which has an outstanding influence on the output of maximum power for wind turbines. The blade of the wind turbine consists of one or more different aerofoil profiles along the blade. The blade can be divided into three main areas: root region, middle region, and tip region. The root region of the blade is built up from aerofoil profiles with a high thickness to chord ratio. Thick aerofoil profiles usually have a lower lift to drag ratio, and therefore, special consideration is required to increase the lift force of the aerofoil for more efficient wind turbine designs (Fuglsang and Bak 2004; Rooij and Timmer 2003). Researches (Anaya-Lara et al. 2011; Schubel and Crossley 2012) stated that wind turbines require a special aerofoil design to increase their power outputs at various flow conditions. In fact, the sudden stall presented on an aerofoil section can also significantly impact on the power output of wind turbines. Therefore, the “gentle stall” is needed for wind turbine aerofoil to avoid such significant losses of energy harvesting (Ahmed 2012).

A great attempt has been done to develop wind turbine technology. National Advisory Committee for Aeronautics (NACA 1949) designed four and five-digit aerofoil series for early modern wind turbines (Hau 2006). These digits classify the geometric profile of the aerofoil as follows: first digit indicates the maximum camber to the chord ratio, second digit refers to the location of the camber in the 10th of the chord, while third and fourth digits together are the percentage maximum thickness to chord ratio (Abbott and Doenhoff 1949). Alternatively, many other energy organizations and universities established different aerofoil profiles for the wind turbine industry such as the Delft University (Timmer and van Rooij 2003), LS, SERI-NREL, and FFA (Burton 2011) and RISO (Fuglsang and Bak 2004).

China has become in the lead throughout the world in the field of wind power generation. The CAS-W1 aerofoil families were designed by the Chinese Academy of Sciences. These aerofoil profiles have thicknesses of 15–25%. (Xu 2010) showed that the CAS-W1-250 aerofoil has a maximum thickness of 25% and a TE thickness of 0.6%. At Re = 3106, the aerofoil has very good aerodynamic characteristics with C_L,max of 1.7 at 15° in clean conditions and 1.66 with LE roughness at the same conditions. At the design angle of attack of 6°, the L/D is equal to 157.6. (Chandrala et al. 2013) selected a single NACA 0018 aerofoil to make a rotor blade. The CFD analysis was carried out using ANSYS CFX software at different angles of attack at a constant wind speed of 32 m/s. The velocity and pressure distributions at various blade angles of attack were determined and they were found to be in good agreement with measurements. The optimal angle of attack at a wind speed of 32 m/s is 10°.

The aerofoil design process is essential in engineering design. It can be categorized into three methods: direct design, inverse design and optimization. In the direct design method (Xu 2010), the specifications of aerofoil design geometry is dictated by the outcome analysis, inferred from wind-tunnel or flight tests. The inverse design method (Vicini and Quagliarella 1997; Petrucci and Filho 2007), determines an aerofoil shape iteratively based on the target surface pressure distribution. There are various optimization methods which are developed over time for aerofoil design such as the conjugate gradient method (XiaoPing et al. 2015), Aquasi-Newton optimization method (Kennelly 1983), the gradient search optimization approach (Zhang et al. 2021), stochastic optimization methods (genetic algorithms (Quagliarella and Cioppa 1995) and simulated annealing (Liu 2005)). A combination between different optimization methods is also available. (Shahrokhi and Jahangirian 2007) applied an optimization approach to the shape of an aerofoil. Two parameterization methods were applied based on the flow characteristics of transonic viscous flow. A Genetic Algorithm was used as the optimization method and the shape of a viscous transonic aerofoil was optimized at high Reynolds number turbulent flow conditions. (Herrmann and Bangga 2019) improved the aerodynamic performance of the thick DU 00-W-401 aerofoil using a genetic algorithm approach for large wind turbine rotors under fixed and free transition conditions. It was found that the aerodynamic performance under fixed conditions was significantly improved, while for the free transition test, the results were slightly improved. (Tesfahunegn et al. 2015) presented a space mapping algorithm to optimize the aerofoil shape considering adjoint sensitivities. The design process showed lower computational costs compared with the high fidelity CFD simulation methods. (Sanaye and Hassanzadeh 2014) carried out a multi-objective optimization approach to maximize the aerofoil lift-to-drag ratio. A genetic algorithm method along with ‘Fuzzy Bellman-Zadeh’ decision-making method was applied as an optimization process. An average increase of 26% in the lift-to-drag ratio of the developed S822 aerofoil was obtained compared to the typical S822 aerofoil.

The main objective of this paper is to design a new aerofoil profile with a high lift-to-drag ratio aerofoil, but also with much softer stall and stable characteristics over various operating conditions that match those experienced by horizontal axis wind turbines (at Reynolds number of 1 × 10⁶). To achieve this, the boundary layer on the upper and lower surfaces of the aerofoil has been investigated, particularly the transitional area and its influence on the aerofoil characteristics. In this paper, the following aspects will be addressed:

1. The relationship between the location of the transitional area and the lift and drag force coefficients for different operating conditions.

2. The impact of the aerofoil geometry (leading edge radius, thickness, and curvature) on the stall characteristics, the lift and drag force coefficients and the transitional area at different angles of attack at Reynolds number of 1 × 10⁶.

3. The influence of the turbulence on the stall behaviour of the new aerofoil characteristics using CFD and measurements.

4. The influence of the combined impact of very high turbulent intensity and surface rough on the performance of the aerofoil at Reynolds number of 1 × 10⁶.

Methodology

The inviscid formulation of XFOIL code

A simple linear-vorticity stream function panel method is used as an inviscid formulation. A two-dimensional inviscid flow field is constructed by the superposition of free stream flow, a vortex sheet of strength γ and a source sheet of strength σ on the aerofoil surface as shown in Figure 1. The trailing edge thickness is modeled by a source panel. The stream function is defined by:

(1)Ψ(x,y)=u∞y−υ∞x+12π∫γ(s) ln r(s;x,y) ds+12π∫σ(s) θ(s;x,y) ds

where r is the magnitude of the vector between the path coordinate S and field points x, y and the vector’s angle θ.

1. Modelling the aerofoil and its wake by panel method along with vorticity and source distributions

2. Details of the trailing edge

Figure 1:

(a) Modelling the aerofoil and its wake by panel method along with vorticity and source distributions, and (b) details of the trailing edge. Figure is adapted from (Drela 1989).

The equations are closed by applying an explicit Kutta condition, in which the source sheet σ is used in the iteration scheme for modelling the displacement effect of the boundary layer. The compressibility effect is taken into account by applying Karman-Tsien compressibility correction for high subsonic conditions.

Computation fluid dynamics (CFD) method

It is important to ensure that the stall characteristics of the developed aerofoil profile are stable and gradual without a sudden drop in the lift force and a sudden increase in the drag force. Therefore, the flow around the new aerofoil was analysed by an open source CFD code, OpenFOAM^® (Greensheilds 2015). In OpenFoam, the standard (k–ω) model is mainly adopted from the work of Wilcox (2008), which uses modifications for the effects of low-Reynolds number, compressibility and shear flow spreading. The standard (k–ω) model is an empirical model which is based on model transport equations for modelling the turbulent kinetic energy (k) and the specific dissipation rate (ω). Both the k and ω equations have been modified over time to improve their accuracy for predicting free shear flows. The general forms of the turbulent kinetic energy (k) and the specific dissipation rate (ω) can be derived from the below transport equations:

(6)∂∂t(ρk)+∂∂xj(ρkui)=∂∂xj(Γk∂k∂xj)+GK+YK+SK

(7)∂∂t(ρω)+∂∂xj(ρωui)=∂∂xj(Γω∂ω∂xj)+Gω+Yω+Sω

where G_k is the generation of turbulent kinetic energy due to mean velocity gradients, G_ω is the generation of ω, Γ_k and Γ_ω are the effective diffusivity of k and ω respectively, Y_k and Y_ω are the dissipation of k and ω due to turbulence, and S_k and S_ω are user-defined source terms.

The domain and mesh are presented in Figures 2 and 3. Higher grid density was made close to the aerofoil wall to ensure modelling the boundary layer with sufficient accuracy. For the computational domain in Figure 2, several grids/meshes have been applied using quadrilateral cells as illustrated in Tables 3 and 4. Each grid was made very dense (with a starting cell height of 1e-4) near the aerofoil surface to capture the pressure gradient at the boundary layer with sufficient accuracy. This is important because the adverse pressure gradient causes the flow separation. In the far-field area of the computational domain, the resolution of the mesh was made a standard grid due to the fact that the flow gradients approach zero away from the aerofoil.

Figure 2:

The domain and mesh used for simulation of the flow past aerofoil.

Figure 3:

Mesh of the computational domain and around NACA 2415 aerofoil by using the structured grid.

Table 3:

Comparisons of the lift and drag force coefficients between experiment (Abbott and Von Doenhoff 2012) and CFD results (NACA2415 aerofoil) for different meshes at 8° angle of attack at Re = 3 × 10⁶.

Grids	Number of cells	Average y+	Experiment		CFD		Error (%)
			C L	C D	C L	C D	ε CL	ε CD
Grid 1	6200	46	1.12	0.012	0.95	0.020	15.17	66.66
Grid 2	10,080	32	1.12	0.012	1.03	0.015	8.03	25
Grid 3	219,600	2.45	1.12	0.012	1.12	0.012	0	0
Grid 4	259,160	0.55	1.12	0.012	1.12	0.012	0	0

Table 4:

Comparisons of the lift and drag force coefficients between experiment (Abbott and Von Doenhoff 2012) and CFD results (NACA2415 aerofoil) for different meshes at 15° angle of attack at attack at Re = 3 × 10⁶.

Grids	Number of cells	Average y+	Experiment		CFD		Error (%)
			C L	C D	C L	C D	ε CL	ε CD
Grid 1	6200	46	1.37	0.021	1.05	0.015	23.35	28.57
Grid 2	10,080	32	1.37	0.021	1.29	0.019	5.83	9.52
Grid 3	219,600	2.45	1.37	0.021	1.34	0.022	2.18	4.76
Grid 4	259,160	0.55	1.37	0.021	1.345	0.020	1.82	4.76

In this approach conservation of mass and momentum equations are discretized according to the standard Gaussian finite-volume integration. These equations are integrated over the volume, whereas the convective and diffusive terms are converted to surface integral with the use of Gauss’ divergence theorem. The volume and surface integrals are approximated by a second-order accurate midpoint scheme. In the latter approach, a linear interpolation is applied to determine the unknown values on the cell faces. Gaussian integration is applied by taking the sum of the values on the cell faces.

Since the turbulent flow field is described by the (k-ω) model, it is necessary to implement a very small time step according to the grid refinement. Therefore, the courant number (C_o = U_∞Δt/Δx), where Δx is the average distance between two cell centroids on the aerofoil wall, and Δt is the time step, was selected to be ≤1 to ensures that fluid particles move from one cell to another within one time step with stable and converged solutions, and optimal error damping properties. The average wall y⁺ distributions over the NACA2415 aerofoil at α = 8^ο and 15^ο at Re = 3 × 10⁶ are also presented in Tables 3 and 4. From Tables 3 and 4, it is obvious that the Grids 1 and 2 give unacceptable error compared to measurements as determined by Eq. (8), while the Grids 3 and 4 give satisfied results. The wall y⁺ distributions for Grids 3 and 4 are illustrated in Figure 4. In this paper, the y⁺ value was chosen to be less than unity for best accuracy as recommended by (Versteeg and Malalasekera 2007). Indeed, the required time step to simulate the flow on the mesh corresponding to Grid 4 is in the order of Δt = O(10⁸) s. Therefore, all simulations were run using Grid 4 with models having a chord length of 1 m at Reynolds number of 1 × 10⁶, which corresponds to wind speed of 15 m/s respectively.

(8)ϵ=|True value–ApproximationTrue value|×100

Figure 4:

y+ distribution over the NACA2415 aerofoil at α = 8^ο at Re = 3 × 10⁶.

Subsonic wind tunnel measurements

The experimental work is carried out in an open-ended subsonic wind tunnel located at the laboratory of the Omar Al-Mukhtar University, Derna, Libya. As seen from Figure 5, the subsonic wind tunnel is provided with a fan driven by a DC variable speed motor unit located downstream of the working section. The DC motor allows regulating the number of revolutions and therefore, the airspeed between 0 and 26 m/s. The air enters the test section through a designed contraction followed by an aluminum honeycomb flow straightener to create uniform flow and to break the large vortexes into small vortexes within the working test section. There are also inclined and multi-tube manometers to measure the air speed and static pressure in the wind tunnel. The dimensions of the wind tunnel are 2.98 m in length, 1.83 m in height, 0.8 m in width, and the test section is 30 × 30 ×45 cm. The turbulent intensity of the wind tunnel is around 0.07% at the maximum wind speed with a mean velocity variation of ±0.2 m/s.

Figure 5:

Open-ended subsonic wind tunnel.

The subsonic wind tunnel is provided with a Pitot static tube (see Figure 6) to determine the air speed in the working test section. It has a 4 mm diameter stainless steel tube and of Prandtl design with a negligible correction up to yaw angles of at least 5°. The Pitot static tubes are connected to the multi-tube manometer by two tappings in order to measure the difference in pressure between the total and static tappings. The air speed in the working section is calculated by applying Bernoulli’s equation as follows:

(9)U∞=2 (ΔP)ρair

where ρ_w is the density of air (kg/m³) and ΔP is the difference in pressure between the total and static tappings (N/m²). ΔP is measured using the tube manometer as follows:

(10)ΔP=ρgΔh

where ρ is manometer fluid density (kg/m³), (water is used), g is gravitational constant (9.81 m/s²), and Δh is the true difference in manometer heights.

Figure 6:

Pitot static tube.

The pressure distribution on the aerofoil can be expressed in terms of a dimensionless lift coefficient C_P which, in turn, can be determined experimentally or numerically.

(11)Cp=pi−p∞12ρU∞2

where p_i is the surface pressure measured at location i on the surface, p_∞ is the pressure in the free stream, ρ is air density, and U_∞ is the free-stream velocity.

One important consideration to take into account in the subsonic wind tunnel is the presence of the solid tunnel boundaries which can cause variation in the flow field from the free-air flow. In order to avoid bounded air stream and inaccurate measurements, a set of corrections caused by the wind tunnel walls must be applied. In the current study, the phenomena of horizontal buoyancy, solid blockage, wake blockage, and stream curvature were all corrected according to the methods of (Allen 1944; Barlow, Rae, and Alan 1999).

Generating turbulence in wind tunnel experiments

Generation of turbulence in wind tunnels is important for studying wind turbine aerodynamics (Sicot et al. 2008). A turbulent inlet is usually generated using three families of grids: passive, active, and fractal. In this paper, the passive grid technique (Taylor 1935) is selected to generate wind tunnel turbulence. It should be noted that the designed turbulent inflow depends mainly on (1)- the width b of the bars, (2)- the mesh size M presenting the distance between the centerline of the bars, and (3)- the downstream distance x to the grid as illustrated in Figure 7. It was shown by (Vickery 1966) that the optimal mesh size is at M = L/8 where L indicates the test section length. There are many empirical relations for turbulence characteristics as presented in Table 5. However, these empirical relationships depend on constants, evaluated in a certain wind tunnel. In this research work, it was found that the most reliable formula for estimating the turbulent intensity is given by (Laneville 1973) as follows:

(12)I=2.54(x/b)−8/9

Figure 7:

Schematic of the grid with bar b and mesh M size.

Table 5:

Empirical relations for turbulence characteristics (Roach 1987).

Empirical expression	I = A(x/b)	I = AB(x/b)	L/b = C(x/b)
Constants	A=1.13	B=0.89	C=0.20

In this paper, it is desired to investigate the stall characteristics at turbulent intensities of 0.50 and 0.75%. Therefore, the passive grid was designed carefully to generate the required turbulent intensities at the inlet of the wind tunnel. Based on the optimal mesh size of M = L/8 and the test section length (L) of 45 cm, M is determined to be 0.0563 m. The width b of the bars, and the downstream distance x to the grid have been varied and the corresponding turbulent intensities are determined by Eq. (12) as seen from Table 6.

Table 6:

Geometry of the designed grids and the relevant turbulent intensities.

Grid	x (m)	M (m)	b (m)	I (%)
1	0.4	0.0563	0.005	5.1
2	0.8	0.0563	0.015	7.56

The passive grids in Table 6 were fabricated and placed at the inlet of the wind tunnel to generate the required turbulent intensities of 0.50 and 0.75% respectively. The wind speed was then measured at Reynolds number of 1 × 10⁶ over a period of 1 h. The measured standard deviations of the wind speed at Re = 1 × 10⁶ are 0.72 for I = 0.50% and 1.1 for I = 0.75%. These measured wind speeds will be used to analyse the impact of the turbulent intensity on the transitional area occurring on the upper surface of the aerofoil at different angles of attack and at Reynolds number of 1 × 10⁶.

The integral length scale (l) is determined by Taylor’s frozen eddy hypothesis as given by:

(13)l=T×U∞

where the integral time scale (T) is calculated by:

(14)T=∫f(τ)dτ

The auto-correlation function of time (f(τ)) is

(15)f(τ)=1u′2∫u(t)u(t−τ)dt

where u is the velocity fluctuation term and u’ is the standard deviation of the velocity fluctuation.

Results and discussion

This section presents the results conducted on the aerofoil profiles developed by the NACA (1949). NACA aerofoil profiles were compared to select the aerofoil profiles with the highest lift-to-drag ratios. From the comparisons, it was found that NACA 2415 has the highest lift-to-drag ratio of 103.702 at Re = 1 × 10⁶ as seen from Table 7, which depicts only results of some selected profiles for illustration of the comparisons. NACA 2415 has a maximum camber 2% of the chord length located at 40% chord from the leading edge with a maximum thickness 15% of the chord at 30% chord from the leading edge as shown in Figure 8. When the suitable aerofoil profile was selected, the baseline aerofoil (NACA 2415) was optimized to further increase the lift-to-drag ratio using the optimization approach described early. The optimized aerofoil profile has the following specifications: a maximum camber 5% of the chord length located at 40% chord from the leading edge with a maximum thickness 15% of the chord at 30% chord from the leading edge as illustrated in Figure 8.

Figure 8:

The geometry of the baseline aerofoil (NACA 2415), and the optimized aerofoil profiles.

Figure 9 illustrates comparisons for the lift force coefficients and lift-to-drag ratios with angles of attack between the baseline aerofoil (NACA 2415) and the optimized aerofoil shape profile at Re = 1 × 10⁶ for clean condition (N_Crit = 9). In these results, the optimized aerofoil shape profile shows higher lift force coefficient and lift-to-drag ratio. With this improvement, the optimized aerofoil shape profile shows an average increase of more than 25% in the lift-to-drag ratio. These results show clearly that the curved surfaces generate higher lift forces than the flat or less curved surfaces. In Figure 9, it can also be seen that the decrease in the lift force coefficient at stall for the optimized aerofoil shape profile is more gradual and smooth compared with that for the baseline aerofoil (NACA 2415). In the attached flow regime, the lift force coefficient of the optimized aerofoil shape profile is higher than that for the baseline aerofoil (NACA 2415). The change from the attached flow regime to the post stall flow regime occurs earlier for the baseline aerofoil (NACA 2415) at static stall angle of attack of 15⁰, while it occurs at 16⁰ for the optimized aerofoil shape profile. This delay is mainly related to changing in the curvature of the aerofoil and the boundary layer flow separation.

Figure 9:

Comparisons between characteristics of the baseline aerofoil (NACA 2415) and the optimized aerofoil shape profiles at Re = 1 × 10⁶ for clean condition (N_Crit = 9).

Figure 10:

Comparisons of the time averaged pressure coefficients between measurements with error bars and CFD results at Re = 1 × 10⁶.

1. Investigation of the stall characteristics: In the previous comparisons, it was depicted that the baseline aerofoil (NACA 2415) is the most suitable aerofoil profile. Therefore, the focus will be only on this profile in the following discussion. Table 8 presents the transitional area, determined by XFOIL code for clean conditions (N_Crit = 9), on the upper surfaces of the baseline aerofoil (NACA 2415) and the optimized aerofoil shape profiles at different angles of attack at Re = 1 × 10⁶. From Table 8, it can be seen that the transitional point on the upper surface of the optimized aerofoil shape profile is delayed further downstream towards the trailing edge by increasing the angle of attack compared to the baseline aerofoil (NACA 2415), except at the deep stall flow regime (α = 20⁰), where the transitional point comes closer to the leading edge of the optimized aerofoil shape profile. Table 8 shows that the optimized aerofoil has a larger laminar boundary layer region than the baseline aerofoil (NACA 2415), which benefits the lift force coefficient.

1. Impact of turbulence on the aerofoil characteristics: The strategy of the generation of turbulence in the wind tunnel was presented before (see Figure 7 and Table 6). A qualitative comparison of the CFD and experimental results of the surface pressure coefficient is illustrated in Figure 10. In these results, the surface pressure distributions at each turbulent intensity at Re = 1 × 10⁶ were measured over a period of 1 h. Each has been repeated 10 times for several days to have enough data points. In Figure 10, the error bars with measurements indicate around two standard deviations. In this comparison, the CFD results are in good agreement with measurements at all angles of attack and turbulent intensities.

Figure 11 shows velocity counters of the CFD results without turbulence for the optimized aerofoil shape profile at Reynolds number of 1 × 10⁶, while Figures 12 and 13 illustrate the same but for turbulent intensities of 0.50 and 0.75%. Form these results, it could be seen that the level of turbulence on the flow separation and the boundary layer momentum thickness has no significant influence in the attached flow regimes (at angles of attack of 5° and 10°) as shown in Figures 11(a, b)–13(a, b). Therefore, it can be concluded that the boundary layer momentum thickness and flow separation are constant and independent on pressure gradient in the attached flow regime for different turbulent intensities. However, different instability mechanisms are observed at near stall (α = 15°) and post stall (α = 20°) flow regimes for different turbulent intensities (see Figures 11(d, c)–13(d, c)). From these results, it is obvious that the structure of the separated shear layer and the transition criteria are slightly unstable for various turbulent intensities.

Figure 11:

Contours of velocity magnitude for the optimized aerofoil shape profile at Re = 1 × 10⁶ (no turbulence).

Figure 12:

Contours of velocity magnitude for the optimized aerofoil shape profile at Re = 1 × 10⁶ (I = 0.50%).

Figure 13:

Contours of velocity magnitude for the optimized aerofoil shape profile at Re = 1 × 10⁶ (I = 0.75%).

To demonstrate the aerodynamic performance of the optimized aerofoil profile, a comparison with the baseline aerofoil (NACA 2415) is conducted at different turbulent intensities and also for smooth and rough surface flows at the same working conditions (Re = 1 × 10⁶, Ma = 0.2) as depicted in Figures 14 and 15. The roughness influence is important to investigate given that the loss of wind turbine power is due to the increased leading edge roughness (caused by insects, sand, salt, and rain over time) on the rotor blades, yielding to reduction in the aerodynamic performance. This impact is taken into account in the XFOIL code by forcing the transition of the boundary layer at 0.01 on the suction side and 0.1 on the pressure side. The maximum lift-to-drag ratio of the optimized aerofoil is 119.2 for the working condition of the smooth surface. This occurs at angle of attack of 6°, while for rough surface condition, the maximum lift-to-drag ratio of the optimized aerofoil is 57.3 at angle of attack of 6°. Although the optimized aerofoil profile is sensitive to the leading edge roughness condition, it still has a good stall characteristic. For the other range of angles of attack, the optimized aerofoil shows higher lift force coefficient and lift-to-drag ratio compared to the baseline aerofoil (NACA 2415) for both free and fixe transitions. It could be observed that the leading edge roughness does not influence on the lift force coefficient slope, but it leads to a reduction in the maximum lift force coefficient and an increase in the minimum drag force coefficient, and subsequently in the lift-to-drag ratio.

Figure 14:

Lift force coefficients and lift-to-drag ratio for clean (free transition) and rough (fixed transition) operating conditions at Reynolds number of 1 × 10⁶.

Figure 15:

Lift force coefficients and lift-to-drag ratio for different turbulent intensities at Reynolds number of 1 × 10⁶-free transition test case.

To analyse the influence of turbulent intensities, the XFOIL code is used with the following specifications: N_Crit = 4.25 for turbulent intensity of 0.50%, and N_Crit = 3.25 for turbulent intensity of 0.75% as shown in Figure 15. As it can be seen, when the level of turbulence is increased from 0.50 to 2.96%, the lift force coefficient is further reduced due to the increase in the turbulent kinetic energy in the flow separation zone of boundary layer and subsequently, the drag force coefficient is increased. At high turbulent intensity of 2.96%, the stall angle of attack occurs earlier than that at low turbulent intensities of 0.50 and 0.75%. This indicates that the transition of the boundary occurs close to the leading edge of the aerofoil under such rough conditions. The effect of the turbulent intensity also causes a small increase in the drag coefficient until the stall angle is reached, after which the drag force coefficient is increased suddenly. This is mainly due to the increase in the pressure drag force as a result of flow separation. In the attached flow regime, it can be observed that the slope of the lift force coefficients is almost the same for all turbulent intensities. When the static stall angle of attack is reached, the lift force coefficient decrease smoothly, indicating smooth and gradual stall characteristics.

An interesting feature, which has not been widely studied in the literature, is to investigate the combined influence of both turbulent intensity and surface roughness on the aerodynamic performance of the aerofoil. In the present work, this was achieved in the XFOIL code by setting a fixed transition to 0.01 and 0.1 on the upper and lower surfaces respectively and N_Crit = 0.01 (I = 2.96%) as illustrated in Figure 16. For such rough working conditions, it is obvious that the optimized aerofoil performs better than the baseline aerofoil (NACA 2415) under the same flow conditions in terms of the lift force coefficient, lift-to-drag ratio, and stall characteristics. From these results, the combined impact of turbulent intensity and surface roughness yields to further drops in the lift force coefficient and the lift-to-drag ratio compared to those shown in Figures 14 and 15.

Figure 16:

Lift force coefficients and lift-to-drag ratio for rough working conditions. In these results the transition is forced at 0.01 and 0.1 on the upper and lower surfaces respectively, while turbulent intensity of 2.96% corresponding to N_Crit = 0.01 is considered.

Conclusions

The current paper was conducted to develop a new aerofoil profile with a high lift-to-drag ratio, smooth and gradual stall for wind turbine applications. The proposed optimization method is based on the use of the experimental Box–Behnken design and XFOIL code, computational fluid dynamics (CFD), and measurements to achieve the purpose. From this work, conclusions are drawn as follows:

1. A new aerofoil profile was developed from the geometric profile of the NACA 2415 aerofoil with the following specifications: a maximum camber of 5% chord located at 40% chord from the leading edge with a maximum thickness of 15% chord located at 30% chord from the leading edge.

2. The newly aerofoil profile provides more than 25% average increase in the lift-to-drag ratio compared to the baseline aerofoil (NACA 2415) under different working conditions. This shows a confidence in utilizing the present methodology for different aerofoils shapes of wind turbine blades.

3. The aerofoil thickness has a dramatic impact on the characteristics of the aerofoil. It was found that thicker aerofoil profiles increase the drag forces and decrease the lift forces, leading to a significant decrease in the lift-to-drag ratio. It is also recommended to consider the locations of the maximum thickness and camber of the aerofoil in the future to further improve the performance of the aerofoil profile using the methodology adopted in this paper.

4. To increase the lift-to-drag ratio, it is suggested that the transitional point on the upper surface of the aerofoil is delayed downstream from the leading edge using a suitable strategy.

5. High turbulent intensity and surface roughness increase the minimum drag force coefficient and reduce both the maximum lift force coefficient and the lift-to-drag ratio. Low turbulent intensity and smooth surface lead to a delay in the stall angle of attack and the flow separation point and thereby to an increase in the lift force coefficient.

6. The newly developed aerofoil has smooth and gradual stall characteristics at different operating conditions. The present study suggests the use of this aerofoil at middle span and blade tip sections of wind turbines since its gives good on-design and off-design aerodynamic performances. The blade root sections of a wind turbine are usually made of thick aerofoil profiles to resist the high stresses. Therefore, this needs to be investigated in the future using the strategy applied in this paper.

7. The main advantage of the proposed optimization method is that satisfactory results are obtained with minimum computational costs. This is a very important feature for the design or optimization of wind turbine blades using a low fidelity accurate flow solver.