1 Introduction
Improving the stall characteristics of aerofoils has been one of the major areas of study in the past. Both active and passive methods, including geometry modifications of trailing and leading edges, use of multi-element aerofoils, suction and blowing mechanisms, Gurney flaps, etc. have been used. More recently, the flippers of a peculiar creature have drawn attention of many researchers. This creature is known as the Humpback whale and is very popular with Australian and American ‘whale watchers’.
Fish & Battle (Reference Fish and Battle1995) performed a very detailed study on the Humpback whale, or Megaptera novaeangliae, that owns the record of having the longest flipper among all cetaceans. Its particularity is the presence of large protuberances or tubercles located at the leading edge of the flipper. It is well known by marine biologists that the Humpback whale also has a unique feeding behaviour which demands high manoeuvrability from the animal. Despite its big dimensions (the whale studied by Fish & Battle (Reference Fish and Battle1995) was 9 m long), the humpback whale is capable of performing high speeds and sharp U-turns. For these manoeuvres the whale orientates its flippers towards a high angle of attack (AoA) while still maintaining the lift. It has been reported in previous literature that, as strakes, the tubercles may not increase the maximum lift of the flippers but they help to maintain a certain level of lift while the whale is turning at a high AoA.
Following the investigations of Fish & Battle (Reference Fish and Battle1995) several researchers continued studying the effects of the undulated or wavy leading-edge (WLE) geometry. Different results were obtained depending on whether a wing tip was involved or an infinite span was used. In the cases with a wing tip, it was routinely found that at high angles of attack the WLEs managed to confine the wing-tip flow separation to the outboard of the wing resulting in an increased the stall angle and maximum lift coefficient (Fish & Battle Reference Fish and Battle1995; Miklosovic et al. Reference Miklosovic, Murray, Howle and Fish2004; Pedro & Kobayashi Reference Pedro and Kobayashi2008; Ozen & Rockwell Reference Ozen and Rockwell2010; Weber et al. Reference Weber, Howle, Murray and Miklosovic2011; Yoon et al. Reference Yoon, Hung, Jung and Kim2011; Guerreiro & Sousa Reference Guerreiro and Sousa2012). Additionally, Guerreiro & Sousa (Reference Guerreiro and Sousa2012) reported that the performance of an undulated wing increased with its aspect ratio. On the other hand, in the cases with an infinite span, the results show that WLEs tend to result in an earlier stall than conventional straight leading edges (SLEs) although the stall process is less abrupt and the post-stall performance (lift-to-drag ratio) may be higher (Johari et al. Reference Johari, Henoch, Custodio and Levshin2007; Miklosovic, Murray & Howle Reference Miklosovic, Murray and Howle2007; Hansen, Kelso & Dally Reference Hansen, Kelso and Dally2011; Favier, Pinelli & Piomelli Reference Favier, Pinelli and Piomelli2012; Zhang, Wang & Xu Reference Zhang, Wang and Xu2013; Rostamzadeh et al. Reference Rostamzadeh, Hansen, Kelso and Dally2014; Skillen et al. Reference Skillen, Revell, Pinelli, Piomelli and Favier2015).
The effect of Reynolds number on the flow over the WLE geometry has also been studied in the past, indicating that the performance of WLEs is less sensitive to Reynolds numbers in comparison to the SLE cases (Guerreiro & Sousa Reference Guerreiro and Sousa2012). Moreover, the improvement of post-stall performance with WLEs was demonstrated over a wide range of Reynolds numbers (Fish & Battle Reference Fish and Battle1995; Miklosovic et al. Reference Miklosovic, Murray, Howle and Fish2004; Johari et al. Reference Johari, Henoch, Custodio and Levshin2007; Miklosovic et al. Reference Miklosovic, Murray and Howle2007; Pedro & Kobayashi Reference Pedro and Kobayashi2008; Hansen et al. Reference Hansen, Kelso and Dally2011; Weber et al. Reference Weber, Howle, Murray and Miklosovic2011; Yoon et al. Reference Yoon, Hung, Jung and Kim2011; Guerreiro & Sousa Reference Guerreiro and Sousa2012; Zhang et al. Reference Zhang, Wang and Xu2013). There were only a few cases where the WLE cases showed an improvement in the pre-stall performance in a finite-wing configuration at a low Reynolds number (Miklosovic et al. Reference Miklosovic, Murray, Howle and Fish2004; Guerreiro & Sousa Reference Guerreiro and Sousa2012; Zhang et al. Reference Zhang, Wang and Xu2013). Nonetheless, it is uniformly agreed that WLEs tend to provide smoother stall characteristics regardless of the Reynolds number or wing configuration used.
One of the common characteristics of undulated wings is to promote a higher pressure at the peaks and a lower pressure at the troughs (Pedro & Kobayashi Reference Pedro and Kobayashi2008; Yoon et al. Reference Yoon, Hung, Jung and Kim2011; Zhang et al. Reference Zhang, Wang and Xu2013; Skillen et al. Reference Skillen, Revell, Pinelli, Piomelli and Favier2015; Hansen et al. Reference Hansen, Rostamzadeh, Kelso and Dally2016). It has been shown that the low pressure at the troughs is the consequence of an accelerated flow channelled between two adjacent peaks. Such pressure variations across the span promotes a secondary cross-flow that results in a stronger flow exchange/mixing. There exists diverse observations with regard to the spanwise pressure distribution at high angles of attack. Skillen et al. (Reference Skillen, Revell, Pinelli, Piomelli and Favier2015) showed a spanwise-periodic pattern where low-pressure pockets of a similar strength were located behind every trough. Other researchers (Custodio Reference Custodio2007; Dropkin et al. Reference Dropkin, Custodio, Henoch and Johari2012; Rostamzadeh et al. Reference Rostamzadeh, Hansen, Kelso and Dally2014) reported that the spanwise flow length scale was no longer equal to the WLE wavelength, i.e. the low-pressure pockets were found to group together in a collocated fashion.
There exists a consensus among Johari et al. (Reference Johari, Henoch, Custodio and Levshin2007), Hansen et al. (Reference Hansen, Kelso and Dally2011), Weber et al. (Reference Weber, Howle, Murray and Miklosovic2011), Favier et al. (Reference Favier, Pinelli and Piomelli2012), Rostamzadeh et al. (Reference Rostamzadeh, Hansen, Kelso and Dally2014) and Skillen et al. (Reference Skillen, Revell, Pinelli, Piomelli and Favier2015) that the flow past a WLE seems to easily separate at the troughs while it remains attached at the peaks. The flow past a trough which starts with a low pressure encounters a large adverse pressure gradient and is therefore prone to an early separation. Studies at low and medium Reynolds numbers has shown that the early flow separation in the trough area is followed by a reattachment further downstream forming a laminar separation bubble (LSB). At high angles of attack, some of the LSBs tend to break down due to the high adverse pressure gradient there, while others seem to group together. These events leave footprints of the low-pressure pockets near the leading edge. Hansen et al. (Reference Hansen, Rostamzadeh, Kelso and Dally2016) studied the origin of LSB in the trough area at a low Reynolds number (
$Re_{\infty }=2230$
) by using particle image velocimetry (PIV) and a steady-state Reynolds-averaged Navier–Stokes (RANS) calculation. Due to the very low Reynolds number, the LSB, which they refer to as the vorticity canopy, covered almost the entire chord length. They compared the LSB zone to an ‘owl face’ of the first kind (Perry & Hornung Reference Perry and Hornung1984).
Another commonly observed feature in the study of WLEs is the counter-rotating streamwise vortices which seem to originate from the trough areas (Johari et al.
Reference Johari, Henoch, Custodio and Levshin2007; Pedro & Kobayashi Reference Pedro and Kobayashi2008; Hansen et al.
Reference Hansen, Kelso and Dally2011; Weber et al.
Reference Weber, Howle, Murray and Miklosovic2011; Yoon et al.
Reference Yoon, Hung, Jung and Kim2011; Favier et al.
Reference Favier, Pinelli and Piomelli2012; Rostamzadeh et al.
Reference Rostamzadeh, Hansen, Kelso and Dally2014; Skillen et al.
Reference Skillen, Revell, Pinelli, Piomelli and Favier2015; Hansen et al.
Reference Hansen, Rostamzadeh, Kelso and Dally2016). Favier et al. (Reference Favier, Pinelli and Piomelli2012) stated that the vortices were created due to the spanwise velocity gradient that triggers Kelvin–Helmholtz instabilities. These instabilities produce vortices in the direction normal to the surface and are immediately tilted and convected by the mean flow. On the other hand, Rostamzadeh et al. (Reference Rostamzadeh, Hansen, Kelso and Dally2014) and Hansen et al. (Reference Hansen, Rostamzadeh, Kelso and Dally2016) suggest that the appearance of the streamwise and vertical vortices is better explained based on the development of Prandtl’s secondary flow of the first type, where the skewness of the flow diverts the initial spanwise vortices in the other directions. Hansen et al. (Reference Hansen, Kelso and Dally2011) reported that the interaction between these vortices was governed by the WLE wavelength
$\unicode[STIX]{x1D706}_{LE}$
. They also reported that, similar to vortex generators, the sweep angle
$h_{LE}/\unicode[STIX]{x1D706}_{LE}$
was the key parameter to define the vortex structure. They compared the WLE amplitude with the vortex generator’s device height.
Reduced unsteady fluctuations in the aerodynamic force has also been observed in recent years (Favier et al. Reference Favier, Pinelli and Piomelli2012; Lau, Haeri & Kim Reference Lau, Haeri and Kim2013; Skillen et al. Reference Skillen, Revell, Pinelli, Piomelli and Favier2015), which makes the WLE geometry very promising for applications operating in highly disrupted flow conditions. In particular, Lau et al. (Reference Lau, Haeri and Kim2013) has discovered that WLEs reduce wall pressure fluctuations when an aerofoil is subjected to a vortical gust in the mean flow. The reduction in unsteady fluctuations is achieved by an enhanced phase interference in the leading-edge response (due to the undulation) and a lower effective impact velocity (normal to the surface) along the hill side of the WLE. The same principles apply to the reduction of aerofoil–turbulence interaction noise as investigated by Kim & Haeri (Reference Kim and Haeri2015), Narayanan et al. (Reference Narayanan, Chaitanya, Haeri, Joseph, Kim and Polacsek2015), Kim, Haeri & Joseph (Reference Kim, Haeri and Joseph2016).
Despite the fact that the topic of aerofoils with a WLE has been covered by several groups of researchers in the past, the understanding as to how WLEs improve post-stall performance is still underdeveloped and in debate. Miklosovic et al. (Reference Miklosovic, Murray and Howle2007), Pedro & Kobayashi (Reference Pedro and Kobayashi2008) and Ozen & Rockwell (Reference Ozen and Rockwell2010) found that WLEs acted as spanwise fences that prevented wing-tip flow separation from spreading inboard. Others speculated that an enhanced momentum exchange in the boundary layer was promoted by the WLE-induced vortical structures (Fish & Battle Reference Fish and Battle1995; Miklosovic et al. Reference Miklosovic, Murray, Howle and Fish2004; Pedro & Kobayashi Reference Pedro and Kobayashi2008; Hansen et al. Reference Hansen, Kelso and Dally2011). Custodio (Reference Custodio2007) and Miklosovic et al. (Reference Miklosovic, Murray and Howle2007) claimed that the counter-rotating streamwise vortices usually laid on the top of the low-pressure zones were reminiscent of the vortex-lift mechanism of delta wings. However, these claims and speculations have not fully been supported by high-fidelity flow data and are also limited in a time-averaged domain. Therefore, this paper aims to provide a more detailed investigation into the post-stall behaviours of WLEs in terms of both time-averaged and unsteady aerodynamics by using high-resolution large-eddy simulations. This paper is structured in the following order. In § 2, the current computational set-up is introduced. In § 3, the time-averaged numerical solutions are analysed to understand the effects of WLEs on post-stall aerodynamic performance. The unsteady effects of WLEs on vortex dynamics and force fluctuations are discussed in § 4. Finally, concluding remarks are given in § 5.
2 Problem description and methodology
The aerofoil section considered in the present study is NACA0021 that has previously been used to investigate the effect of undulated leading edges by Hansen et al. (Reference Hansen, Kelso and Dally2011), Rostamzadeh et al. (Reference Rostamzadeh, Hansen, Kelso and Dally2014) and Skillen et al. (Reference Skillen, Revell, Pinelli, Piomelli and Favier2015). The chord-based Reynolds number and the incidence angle (angle of attack or AoA) are set to
$Re_{\infty }=1.2\times 10^{5}$
and
$\unicode[STIX]{x1D6FC}=20^{\circ }$
, respectively, as previously studied in the literature. The free-stream Mach number is set to
$M_{\infty }=0.3$
in this work, at which compressibility effects are insignificant and the computational cost is fairly comparable to that of an incompressible solution. The current computational study employs a high-resolution implicit large-eddy simulation (ILES) method based on a wavenumber-optimised discrete filter. The filter is applied directly to the solution (conservative variables) at every time step and acts as an implicit subgrid scale (SGS) model that enforces dissipation of scales smaller than the filter cutoff wavelength. Garmann, Visbal & Orkwis (Reference Garmann, Visbal and Orkwis2013) performed an extensive analysis of the ILES technique in comparison to more conventional ones with an explicit SGS model and concluded that, subject to an appropriate grid resolution, ILES simulations are capable of correctly capturing the flow physics of a flow such as the one being studied in the current study.
2.1 Governing equations and numerical methods
In this work, the full compressible Navier–Stokes equations (non-dimensionalised) are used, which can be expressed in a conservative form transformed onto a generalised coordinate system as
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\left(\frac{\boldsymbol{Q}}{J}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}_{i}}\left(\frac{\boldsymbol{E}_{j}}{J}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}_{i}}{\unicode[STIX]{x2202}x_{j}}\right)=\frac{M_{\infty }}{\mathit{Re}_{\infty }}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}_{i}}\left(\frac{\boldsymbol{F}_{j}}{J}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}_{i}}{\unicode[STIX]{x2202}x_{j}}\right)-\frac{\boldsymbol{S}}{J},\end{eqnarray}$$
for
$i=1$
, 2, 3; and,
$j=1$
, 2, 3. The Reynolds number is
$Re_{\infty }=\unicode[STIX]{x1D70C}_{\infty }U_{\infty }L_{c}/\unicode[STIX]{x1D707}_{\infty }$
based on the mean chord length (
$L_{c}$
) and the free-stream velocity (
$U_{\infty }=\sqrt{u_{\infty }^{2}+v_{\infty }^{2}+w_{\infty }^{2}}$
); and, the Mach number is
$M_{\infty }=U_{\infty }/a_{\infty }$
where
$a_{\infty }=\sqrt{\unicode[STIX]{x1D6FE}p_{\infty }/\unicode[STIX]{x1D70C}_{\infty }}$
is the ambient speed of sound. The constituting terms in (2.1) are given as follows:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{Q}=[\unicode[STIX]{x1D70C},\unicode[STIX]{x1D70C}u,\unicode[STIX]{x1D70C}v,\unicode[STIX]{x1D70C}w,\unicode[STIX]{x1D70C}e_{t}]^{\text{T}},\\ \boldsymbol{E}_{j}=[\unicode[STIX]{x1D70C}u_{j},(\unicode[STIX]{x1D70C}uu_{j}+\unicode[STIX]{x1D6FF}_{1j}p),(\unicode[STIX]{x1D70C}vu_{j}+\unicode[STIX]{x1D6FF}_{2j}p),(\unicode[STIX]{x1D70C}wu_{j}+\unicode[STIX]{x1D6FF}_{3j}p),(\unicode[STIX]{x1D70C}e_{t}+p)u_{j}]^{\text{T}},\\ \boldsymbol{F}_{j}=[0,\unicode[STIX]{x1D70F}_{1j},\unicode[STIX]{x1D70F}_{2j},\unicode[STIX]{x1D70F}_{3j},u_{i}\unicode[STIX]{x1D70F}_{ji}+q_{j}]^{\text{T}},\\ \unicode[STIX]{x1D70F}_{ij}=\unicode[STIX]{x1D707}\left({\displaystyle \frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}}+{\displaystyle \frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}}-{\displaystyle \frac{2}{3}}\unicode[STIX]{x1D6FF}_{ij}{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{i}}}\right),\quad q_{j}={\displaystyle \frac{\unicode[STIX]{x1D707}}{(\unicode[STIX]{x1D6FE}-1)Pr}}{\displaystyle \frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}x_{j}}},\end{array}\right\}\end{eqnarray}$$
where
$x_{j}=\{x,y,z\}$
and
$\unicode[STIX]{x1D709}_{i}=\{\unicode[STIX]{x1D709},\unicode[STIX]{x1D702},\unicode[STIX]{x1D701}\}$
stand for Cartesian and generalised coordinates, respectively. The Jacobian determinant of the coordinate transformation (from Cartesian to body-fitted) is given by
$J^{-1}=|\unicode[STIX]{x2202}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702},\unicode[STIX]{x1D701})/\unicode[STIX]{x2202}(x,y,z)|$
(Kim & Morris Reference Kim and Morris2002). The vector
$\boldsymbol{Q}$
contains the conservative variables where
$u_{j}=\{u,v,w\}$
are the velocity components and
$e_{t}=p/\left[(\unicode[STIX]{x1D6FE}-1)\unicode[STIX]{x1D70C}\right]+u_{j}u_{j}/2$
is the total energy. In this work,
$\unicode[STIX]{x1D6FE}=1.4$
and the Prandtl number is
$\mathit{Pr}=0.71$
for air. The convective, viscous and heat fluxes are represented by
$\boldsymbol{E}_{j}$
,
$\boldsymbol{F}_{j}$
and
$q_{j}$
, respectively, and
$\unicode[STIX]{x1D70F}_{ij}$
is the stress tensor. The extra source term
$\boldsymbol{S}$
on the right-hand side of (2.1) is non-zero within a prescribed ‘sponge layer’ only (Kim, Lau & Sandham Reference Kim, Lau and Sandham2010; Kim & Haeri Reference Kim and Haeri2015) and is intended to provide a perfectly anechoic condition suppressing numerical reflections from the outer boundaries.
All quantities presented in this paper are dimensionless. All length scales are normalised by
$L_{c}$
, time scales by
$L_{c}/a_{\infty }$
, velocities by
$a_{\infty }$
and pressures by
$\unicode[STIX]{x1D70C}_{\infty }a_{\infty }^{2}$
. In addition, the density, temperature and viscosity are normalised by their own ambient values
$\unicode[STIX]{x1D70C}_{\infty }$
,
$T_{\infty }$
and
$\unicode[STIX]{x1D707}_{\infty }$
, respectively, as usual.
In the current set-up,
$u_{\infty }=U_{\infty }\cos \unicode[STIX]{x1D6FC}$
,
$v_{\infty }=U_{\infty }\sin \unicode[STIX]{x1D6FC}$
and
$w_{\infty }=0$
with the aerofoil’s mean chord placed horizontally on the
$x$
-axis. The flow is gradually accelerated from zero to the desired speed (
$M_{\infty }=0.3$
) by using a moving frame technique for a certain duration of time. This is intended to avoid unwanted ‘added mass’ effects affecting the result and to minimise the initial transition phase of the simulation. The flow acceleration is implemented during the first 30 time units by using a smooth ramping function:
$\boldsymbol{u}_{\mathit{ramp}}=\boldsymbol{u}_{\infty }\sin ^{2}[(\unicode[STIX]{x03C0}/2)\min (t/t_{0},1)]$
where
$t_{0}=30$
.
In this work, the governing equations given above are solved by using high-order accurate numerical methods on a structured multi-block grid system. The flux derivatives in space are calculated based on fourth-order pentadiagonal compact finite difference schemes with seven-point stencils (Kim Reference Kim2007). Explicit time advancing of the numerical solution is carried out by using the classical fourth-order Runge–Kutta scheme with a Courant–Friedrichs–Lewy (CFL) number of 0.95. Numerical stability is maintained by implementing sixth-order pentadiagonal compact filters for which the cutoff wavenumber (normalised by the grid spacing) is set to
$0.85\unicode[STIX]{x03C0}$
(Kim Reference Kim2010). In addition to the sponge layers used, characteristics-based non-reflecting boundary conditions (Kim & Lee Reference Kim and Lee2000) are applied at the far boundaries in order to prevent any outgoing waves from returning to the computational domain. Periodic conditions are used across the spanwise boundary planes. No-slip wall boundary conditions are implemented on the aerofoil surface (Kim & Lee Reference Kim and Lee2004). The computation is fully parallelised via domain decomposition and message passing interface (MPI) approaches, where the distributed solution of the pentadiagonal matrix systems is achieved based on Kim (Reference Kim2013).
2.2 Geometry and discretisation of the problem
The NACA0021 aerofoil geometry with a WLE (and a straight sharp trailing edge) and the computational domain used in this work are shown in figure 1. The positions of the leading and trailing edges are given by
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}x_{LE}=-{\displaystyle \frac{1}{2}}+h_{LE}\sin \left({\displaystyle \frac{2\unicode[STIX]{x03C0}z}{\unicode[STIX]{x1D706}_{LE}}}\right)\quad \text{for }-{\displaystyle \frac{L_{z}}{2}}\leqslant z\leqslant {\displaystyle \frac{L_{z}}{2}},\\ x_{TE}={\displaystyle \frac{1}{2}},\quad y_{LE}=y_{TE}=0,\end{array}\right\}\end{eqnarray}$$
where
$h_{LE}$
and
$\unicode[STIX]{x1D706}_{LE}$
are the amplitude and wavelength of the WLE, respectively. The spanwise domain size is denoted by
$L_{z}$
and it is set to a multiple of
$\unicode[STIX]{x1D706}_{LE}$
(from two to eight). In this work,
$h_{LE}=0.015$
and
$\unicode[STIX]{x1D706}_{LE}=0.11$
are chosen from the previous literature. The origin of the coordinate system is located at the centre of the aerofoil’s centre plane. The domain of investigation comprises
$x\in [-7,14]$
in the streamwise direction,
$y\in [-7,14]$
in the vertical direction and
$z\in [-L_{z}/2,L_{z}/2]$
in the spanwise direction. As it can be seen in figure 1 the aerofoil is located closer to the bottom left boundary so the there is enough space to cover the wake flow behind the aerofoil.
Figure 1. The aerofoil geometry and the computational domain used in the current simulations: (a) the side view of the NACA0021 aerofoil with a WLE, (b) the top view with the surface meshes (every fifth point plotted in both directions) and (c) an
$xy$
-plane cross-section of the computational domain showing a contour plot of the mesh aspect ratio (black lines indicating
$AR_{mesh}=2$
).
Based upon a grid dependency test performed (see figure 2), the selected grid (G1 in table 1) consists of
$38.4\times 10^{3}$
cells on each
$xy$
-plane with 800 cells in the streamwise direction and 480 cells in the vertical direction. For every wavelength of the WLE (
$\unicode[STIX]{x1D706}_{LE}$
) 50 cells are used as standard. However, it is reduced to 25 cells per
$\unicode[STIX]{x1D706}_{LE}$
for the largest cases with
$L_{z}=8\unicode[STIX]{x1D706}_{LE}$
(SLE-8 and WLE-8). The total cell count for the largest simulation is
$76.8\times 10^{6}$
. The simulations are run for up to 200 time units (
$t=200$
), which takes around 18 wall-clock hours (in the largest case including data input and output) using 7680 processor cores in the national supercomputer ARCHER.
Figure 2. The time- and spanwise-averaged profiles of pressure and skin-friction coefficient for the WLE-4 case obtained from four different levels of grid resolution listed in table 1.
Figure 3. Distributions of the current surface mesh sizes in wall units (averaged in span).
Table 1. Four different levels of grid resolution (G1–G4) used for a grid dependency test based on the WLE-4 case. The list shows the number of grid cells used in the streamwise (
$n_{\unicode[STIX]{x1D709}}$
), vertical (
$n_{\unicode[STIX]{x1D702}}$
) and spanwise (
$n_{\unicode[STIX]{x1D701}}$
) directions; and the resulting lift and drag coefficients (time averaged).
Since no wall modelling is considered here, a sufficiently high level of near-wall grid resolution is maintained in order to properly resolve the boundary layers throughout the surface. Figure 3 shows the distribution of the current near-wall grid resolution, which is even suitable for the direct numerical simulation (DNS) criteria suggested by Georgiadis, Rizzetta & Fureby (Reference Georgiadis, Rizzetta and Fureby2010):
$10\leqslant \unicode[STIX]{x0394}x^{+}\leqslant 20$
,
$\unicode[STIX]{x0394}y^{+}<1$
and
$5\leqslant \unicode[STIX]{x0394}z^{+}\leqslant 10$
. The relaxed spanwise resolution for the SLE-8 and WLE-8 cases (mentioned above) is still well within the large-eddy simulation (LES) criteria according to Georgiadis et al. (Reference Georgiadis, Rizzetta and Fureby2010). While the grid is stretched outwards from the surface, the near-wake region still contains
$14.4\times 10^{3}$
cells per
$xy$
-plane so that the vortices shed from the aerofoil are well captured. An additional effort has been made in order to keep the mesh aspect ratio as close to unity as possible within the wake region, with all of the wake laying on cells with
$AR_{mesh}<2$
as seen in figure 1.
The simulations are performed with several different sizes of span in order to show and compare their effects on the resulting flow field. A total of seven different cases are presented in this paper and they are labelled herein as SLE-2, SLE-4, SLE-8, WLE-2, WLE-3, WLE-4 and WLE-8. SLE and WLE represent the straight and wavy leading edges, respectively, and the integer numbers indicate how many wavelengths of the unit WLE (
$\unicode[STIX]{x1D706}_{LE}$
) are covered across the span in each case.
3 Time-averaged characteristics
Table 2 shows the time-averaged force coefficients for each simulation presented here compared with available experimental and computational data in the literature. The table covers various cases with different aspect ratios (spanwise domain lengths) of the aerofoil used for the same flow condition,
$Re_{\infty }=1.2\times 10^{5}$
and
$\unicode[STIX]{x1D6FC}=20^{\circ }$
. The present results are compared with the previous measurement by Hansen et al. (Reference Hansen, Kelso and Dally2011) who used 66 WLE wavelengths in span (
$AR=7.26$
). The dataset reveals that there is a clear trend of convergence in the simulation data towards the experimental ones as the aspect ratio increases. The convergence trend and the level of agreement is consistent for
$\langle C_{L}\rangle$
,
$\langle C_{D}\rangle$
and
$\langle C_{L}\rangle /\langle C_{D}\rangle$
alike. Although some uncertainty in the experimental data is implied in figure 4(a) (compared with data from the 1930s), the authors have to support the new data by Hansen et al. (Reference Hansen, Kelso and Dally2011) because the simulation results display a clear convergence towards them in both SLE and WLE cases (shown in figure 4
b).
Table 2. Time-averaged lift and drag coefficients; aerodynamic efficiency (
$\langle C_{L}\rangle /\langle C_{D}\rangle$
); standard deviations of the lift and drag coefficients (
$\unicode[STIX]{x1D70E}_{C_{L}}$
and
$\unicode[STIX]{x1D70E}_{C_{D}}$
); and relative differences of the lift and drag coefficients with respect to the experimental data by Hansen et al. (Reference Hansen, Kelso and Dally2011) (
$\unicode[STIX]{x1D716}_{C_{L}}$
and
$\unicode[STIX]{x1D716}_{C_{D}}$
). All values listed above are obtained at
$Re_{\infty }=1.2\times 10^{5}$
and
$\unicode[STIX]{x1D6FC}=20^{\circ }$
.
Figure 4. Comparisons of time-averaged lift coefficient at
$\unicode[STIX]{x1D6FC}=20^{\circ }$
: (a)
$\langle C_{L}\rangle$
varying with Reynolds numbers and (b) the relative difference (%) of
$\langle C_{L}\rangle$
against Hansen et al. (Reference Hansen, Kelso and Dally2011) at
$Re_{\infty }=1.2\times 10^{5}$
varying with the spanwise domain sizes (aspect ratios) used. The initials represent data from Stack (Reference Stack1931), Jacobs (Reference Jacobs1932), Hansen et al. (Reference Hansen, Kelso and Dally2011) and Skillen et al. (Reference Skillen, Revell, Pinelli, Piomelli and Favier2015).
Comparing the SLE-8 and WLE-8 cases, the current simulation confirms that the WLE geometry produces a higher
$\langle C_{L}\rangle$
and a lower
$\langle C_{D}\rangle$
than the SLE counterpart as reported in the literature. An additional aspect revealed in this work is that unsteady fluctuations in both
$C_{L}$
and
$C_{D}$
are reduced in the WLE cases (rather consistently regardless of the aspect ratio), which supports the previous observations by Favier et al. (Reference Favier, Pinelli and Piomelli2012) and Skillen et al. (Reference Skillen, Revell, Pinelli, Piomelli and Favier2015). The rest of the paper is then devoted to making detailed investigations into the flow development and behaviour on and around the aerofoil in order to understand and explain how the WLE geometry creates the aerodynamic benefits compared to the SLE counterpart.
3.1 Overview of time-averaged flow field
Figure 5 shows that the flow field is highly three-dimensional in the WLE cases, particularly near the leading edge, and confirms that the solution is significantly dependent on the spanwise domain size used. One of the most distinctive features made by the WLE compared to the SLE case is the low-pressure spots behind the troughs which are identified as laminar separation bubbles (LSBs) – to be shown in figure 9. It is manifested that the LSBs tend to group together creating a large zone of low pressure. It appears that the LSB grouping intensifies the low-pressure spots (figure 6 d) and generates a high lift force from there. On the other hand, in each WLE case, there is only one trough area left without an LSB as far as the current simulations are concerned. The isolated trough areas are submerged in a fully separated flow and make a loss in lift force there. In addition, figure 6 indicates that the gap between the lowest and the largest pressure on the suction side becomes larger as the aspect ratio increases, which results in a strong pressure gradient across the span.
Figure 5. Contour plots of time-averaged wall pressure coefficient on the suction side obtained from the WLE-2, WLE-3, WLE-4 and WLE-8 cases compared to the SLE-8 case at
$Re_{\infty }=1.2\times 10^{5}$
and
$\unicode[STIX]{x1D6FC}=20^{\circ }$
. P1 and T1 denote the first (far bottom) peak and trough, respectively.
Figure 6. Profiles of time-averaged wall pressure coefficient along the chord on the suction side obtained from various
$xy$
-plane cross-sections in the WLE cases compared with the SLE-8 case (averaged in span). All cases are at
$Re_{\infty }=1.2\times 10^{5}$
and
$\unicode[STIX]{x1D6FC}=20^{\circ }$
. The locations of the cross-sections (T1, T2, etc.) are indicated in figure 5.
The distribution of wall shear stress of the WLE cases are shown in figure 7 compared with the SLE case. Typical attached flows behind the peaks and the locally separated flows (LSBs) behind the troughs are displayed, as also reported in the literature. More importantly, it is revealed that the grouping of LSBs is followed by a large zone of attached flow. The size of the attached flow zone becomes so substantial in the WLE-8 case that it stretches all the way down to the point close to the trailing edge. Also, it seems that the central LSB becomes dominant and advances downstream as the size of the LSB group grows. Figure 8 shows the profiles of the skin-friction coefficient along the chord at various spanwise locations. It is clearly manifested in the WLE-8 case that the skin friction is reduced inside the LSBs (where the flow direction is reversed) and increased in the reattached flow region in comparison with the SLE case. Meanwhile, the profile at the isolated trough T8 (with a fully separated flow) is almost same as that of the SLE case. Based on the surface data, the skin-friction drag is calculated and compared between the SLE and WLE cases in table 3. The skin-friction drag is also compared with the pressure drag. The table shows that the skin-friction drag is increased in the WLE case which is due to the contribution of the attached flow behind the LSBs. However, it is obvious that the skin-friction drag is only a tiny fraction of the pressure drag.
Figure 7. Contour plots of streamwise wall shear stress on the suction side obtained from the WLE-2, WLE-3, WLE-4 and WLE-8 cases compared to the SLE-8 case at
$Re_{\infty }=1.2\times 10^{5}$
and
$\unicode[STIX]{x1D6FC}=20^{\circ }$
. P1 and T1 denote the first (far bottom) peak and trough, respectively. The locations of
$\langle \unicode[STIX]{x1D70F}_{w_{x}}\rangle =0$
are highlighted by thick black lines.
Figure 8. Profiles of time-averaged skin-friction coefficient along the chord on the suction side obtained from various
$xy$
-plane cross-sections in the WLE cases compared to the SLE-8 case (averaged in span). All cases are at
$Re_{\infty }=1.2\times 10^{5}$
and
$\unicode[STIX]{x1D6FC}=20^{\circ }$
. The locations of the cross-sections (T1, T2, etc.) are indicated in figure 7.
Table 3. Breakdown of the total drag into the pressure (
$C_{Dp}$
) and skin-friction (
$C_{Df}$
) components: time-averaged and fluctuating (standard deviation) parts.
A clear picture of the separation/reattachment points is found in figure 9 by plotting of iso-surfaces of
$\langle u\rangle =0$
(zero time-averaged streamwise velocity). The iso-surfaces are the location where the direction of the flow reverses and hence show interfaces between the attached and separated flow regions. In figure 9 it is visualised that the WLE cases produce a significantly large zone of attached flow downstream of the LSB group. In the SLE case, the leading-edge flow separation without reattachment is most likely to be due to the LSBs broken down in the early stage of flow development.
Figure 9. Iso-surfaces of
$\langle u\rangle =0$
on the suction side obtained from the WLE-2, WLE-3, WLE-4 and WLE-8 cases compared to the SLE-8 case at
$Re_{\infty }=1.2\times 10^{5}$
and
$\unicode[STIX]{x1D6FC}=20^{\circ }$
. The light grey areas indicate attached flow regions and the dark grey surfaces represent separated shear layers.
Figure 10 provides a more detailed look at the near-wall flow structures obtained by applying the linear integral convolution (LIC) technique (Cabral & Leedom Reference Cabral and Leedom1993) to the wall shear-stress vector field. In essence, the result of implementing LIC is similar to an oil-based surface flow visualisation commonly used in wind-tunnel experiments. It shows that the flow from the peaks bends towards the LSBs (low-pressure spots). Some of the flow reverses its direction to enter the LSBs through the side. The rest is joined and pushed downstream (re-energised) by the stream sliding down from the top of the LSB canopy (to be shown in figure 17), which manages to redevelop an attached flow (highly turbulent boundary layer) downstream of the LSBs. Figure 10 also shows some critical points identified according to Perry & Chong (Reference Perry and Chong1987). A common feature is found at the front parts of the LSBs that the internal flow structure is characterised by two saddle points and two unstable foci. Further downstream, the LSBs create more saddle points where the flow reverses its direction. Also, stable nodes are found in the central LSB and a stable focus is identified in each of the skewed satellite LSBs. Unstable foci driving the flow in the interior of the stall cell region can be found as well.
Figure 10. Near-wall flow structure for the WLE-4 case represented by a linear integral convolution of the wall shear stress coloured by the streamwise wall shear stress where the direction of the flow is indicated by the black arrows and several critical points (Perry & Chong Reference Perry and Chong1987) have been identified.
Figure 11 displays streamlines going through a stall cell associated with the fully separated flow. It is viewed from large-scale perspectives that the stall cell acts as if it interconnects the streamlines leaving the aerofoil into the wake region and the others coming back to the aerofoil. This process results in a closed loop of streamlines consisting of the wake recirculation zone. It is envisaged that the size, location and strength of the stall cell have a close relationship with the rear wake, which deserves further investigations particularly in the study of undulated aerofoils since the rear wake is strongly related with the aerodynamic performance.
Figure 11. Streamlines obtained from time-averaged flow data (WLE-2) visualising the flow travelling through the stall cell. The streamlines are colour coded to distinguish those moving away from the aerofoil (red) and the others coming back to it (blue) where the arrow heads indicate the direction of the streamlines. Additional streamlines (grey) are provided on a fixed cross-section to show the flow recirculation in the wake region.
3.2 Origin of streamwise vorticity at the leading edge
It is shown earlier that the flow passing through a trough area either undergoes a full separation immediately or forms an LSB followed by a reattachment downstream. Regardless, a shear layer appears after each trough where the flow separation occurs. The shear layer then rapidly rolls up and produces vortical structures due to the Kelvin–Helmholtz instability. Figure 12 displays some of the vortical structures captured instantaneously comparing the SLE and WLE cases. It can be seen that spanwise vortices are predominant in the SLE case because of the spanwise uniformity in the flow separation at the leading edge. On the other hand, the WLE case exhibits more diverse vortical structures around the LSB due to a high three-dimensionality in the flow there. In particular, a pair of counter-rotating streamwise vortices surrounding the rear part of the LSB is observed. This observation is in line with those reported in the literature.
Figure 12. Instantaneous iso-surfaces of
$Q$
-criterion (the second invariant of the velocity gradient tensor) coloured by streamwise vorticity for SLE-2 and WLE-2 cases: (a)
$Q=1500$
and (b)
$Q=150$
.
As similarly attempted by Hansen et al. (Reference Hansen, Rostamzadeh, Kelso and Dally2016), the origin and development of these streamwise vortical structures may be explained by using the stretching–turning term of the vorticity transport equation:
$$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D714}_{i}}{\text{d}t}=\underbrace{\unicode[STIX]{x1D714}_{j}\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}}_{\substack{ \mathit{Stretching} \\ \mathit{and~turning}}}-\underbrace{\unicode[STIX]{x1D714}_{i}\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{j}}}_{\substack{ \mathit{Compressibility} \\ \mathit{stretching}}}+\underbrace{e_{ijk}\frac{1}{\unicode[STIX]{x1D70C}^{2}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{j}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{k}}}_{\mathit{Baroclinic}}+\underbrace{e_{ijk}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left(\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{km}}{\unicode[STIX]{x2202}x_{m}}\right)}_{\mathit{Viscous~dissipation}}.\end{eqnarray}$$
According to the equation, the time evolution of streamwise vorticity (
$\unicode[STIX]{x1D714}_{x}$
) is related with the turning of the spanwise (
$\unicode[STIX]{x1D714}_{z}$
) and vertical (
$\unicode[STIX]{x1D714}_{y}$
) vortices driven by the gradients of streamwise velocity, i.e.
$\unicode[STIX]{x1D714}_{z}(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z)$
and
$\unicode[STIX]{x1D714}_{y}(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y)$
. Figure 13 briefly describes the turning mechanism. Figure 14 shows that
$\langle \unicode[STIX]{x1D714}_{z}\rangle \langle \unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z\rangle +\langle \unicode[STIX]{x1D714}_{y}\rangle \langle \unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y\rangle$
produces an almost identical picture to that of
$\langle \unicode[STIX]{x1D714}_{x}\rangle$
. This explains the creation of the streamwise vorticity at the leading edge. Hansen et al. (Reference Hansen, Rostamzadeh, Kelso and Dally2016) have provided a similar explanation at a much lower Reynolds number (
$Re_{\infty }=2230$
). The current study unfolds a very different evolution of the leading-edge vortices afterwards (to follow in § 3.3).
Figure 14. Iso-surfaces of time-averaged vorticity magnitude
$\langle |\unicode[STIX]{x1D74E}|\rangle =45$
coloured by (a,e,i) streamwise vorticity; (b,f,j) streamwise turning of spanwise vorticity; (c,g,k) streamwise turning of vertical vorticity; and, (d,h,l) the sum of the turning components. For three different WLEs: (a–d) WLE-2; (e–h) WLE-3; and, (i–l) WLE-4 cases. SSL represents the fully separated shear layer from the isolated trough without an LSB.
Figure 15. Contour plots of time-averaged streamwise vorticity (
$\langle \unicode[STIX]{x1D714}_{x}\rangle$
) for the WLE-4 case at various streamwise locations. The frames cover from the leading edge to the rear of the LSB group.
Figure 14 also shows that
$\langle \unicode[STIX]{x1D714}_{z}\rangle \langle \unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z\rangle$
(streamwise turning of spanwise vorticity) and
$\langle \unicode[STIX]{x1D714}_{y}\rangle \langle \unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y\rangle$
(streamwise turning of vertical vorticity) appear significant around the LSBs when viewed individually but they seem to disappear when put together. It may be due to the following. The sign of
$\langle \unicode[STIX]{x1D714}_{z}\rangle$
is constantly negative (rolled up) where
$\langle \unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z\rangle$
is negative at the left-hand side of an LSB and positive at the right-hand side viewed from the front (see figure 13). In the meantime,
$\langle \unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y\rangle$
is constantly positive on the LSB where
$\langle \unicode[STIX]{x1D714}_{y}\rangle$
is negative at the left-hand side of the LSB and positive at the right-hand side. Therefore the turning terms yield opposite signs that cancel each other around the LSBs. However, the contribution of
$\langle \unicode[STIX]{x1D714}_{z}\rangle \langle \unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z\rangle$
created right at the leading edge remains strong and convects downstream surrounding the LSBs (figure 14
b,f,j).
3.3 Evolution of streamwise vortices behind the leading edge
In order to investigate the evolution of the streamwise vorticity from the leading edge to the rear of the LSBs, contour plots of
$\langle \unicode[STIX]{x1D714}_{x}\rangle$
are provided in figure 15 at various streamwise locations. Up to the position
$x=-0.47$
the thin streamwise vortex (SV) sheets are observed to remain almost undisturbed. Non-uniform features start to appear at
$x=-0.46$
. The SV sheets are detached straight away from the wall and diffused fast at T4 where the trough undergoes a fully separated shear layer (SSL) as mentioned earlier. At the other troughs (T1 to T3), the SV sheets remain close to the wall but they are slightly lifted up from the wall due to the presence of the LSBs (
$x=-0.45$
). Underneath of the uplifted SV sheets, a thin sublayer of streamwise vorticity (rotating in the opposite direction) starts to emerge (
$x=-0.44$
), which takes place at the skewed satellite LSBs (T1 and T3) but not at the central LSB (T2). The footprints of the sublayers mainly driven by
$\langle \unicode[STIX]{x2202}w/\unicode[STIX]{x2202}y\rangle$
are displayed in figure 16 where the opposite sign of the sublayers to that of the initial SV sheets is revealed. Interestingly, the sublayers appear to interact with the uplifted SV sheets to create a strong vortex roll-up at T1 and T3 (
$x=-0.44$
). The core of the rolled-up vortices is connected with the uplifted SV sheets with the same direction of rotation. At
$x=-0.43$
, the uplifted SV sheets are now reaching deep into the adjacent SSLs. This implies that the rolled-up vortex core is connected with the adjacent SSL from where a feeding of streamwise vorticity (as well as turbulent kinetic energy – to follow) is available. The rolled-up vortices remains as the major flow structure but they eventually fade away downstream since the feeding from the SSLs is no longer available.
Figure 16. Contour plots of velocity gradients contributing to the streamwise vorticity at the wall (on the suction side) for the WLE-4 case.
Figure 17. Time-averaged streamlines on the suction side of the WLE-4 case, showing (a) a pair of counter-rotating streamwise vortices surrounding the LSB group (coloured by
$\langle \unicode[STIX]{x1D714}_{x}\rangle$
) and the internal structure of the central LSB (coloured by
$\langle u\rangle$
) in a perspective view; (b) the side view of the central LSB; and, (c) the top view of the central LSB.
It is observed in the current simulations that only one pair of rolled-up (counter-rotating) SVs survives at the furthest troughs from the central LSB. It is presumably because they need feeding from the adjacent SSLs to maintain the streamwise vorticity. It is shown in figure 17(a) that the pair of rolled-up SVs stretches downstream wrapping around the LSB group and effectively forms an interface/buffer preventing the SSLs from penetrating into the LSB group. The two counter-rotating SVs get closer to each other moving downstream before they end up lifting each other up due to the induced upwash velocity between them. The upwash velocity seems to cease the life of the attached flow as well. The structure of the central LSB protected by the rolled-up SV pair is visualised in figure 17(b,c) in detail. It has two (primary and secondary) recirculation zones in it where the primary structure contains a strong circulation to re-energise the flow and create a large zone of attached flow behind it.
An additional investigation is made here with regard to the distribution of turbulent kinetic energy in order to see how the vortical structures evolve with the flow turbulence. The turbulent kinetic energy (TKE) is defined by
$k=(\langle u^{\prime 2}\rangle +\langle v^{\prime 2}\rangle +\langle w^{\prime 2}\rangle )/2$
where
$\boldsymbol{u}^{\prime }=\boldsymbol{u}-\langle \boldsymbol{u}\rangle$
. Figure 18 shows the distribution of TKE for the WLE-4 case at various streamwise locations. At
$x=-0.43$
the plot shows a clear picture of the major flow characteristics, i.e. the central LSB (T2), the skewed satellite LSBs (T1 and T3), the SSL (T4) and the pair of rolled-up SVs at the side edges of the LSB group. The canopies of the LSBs where a thin turbulent shear layer takes place are well captured. A high level of turbulence is manifested in the thick SSL. Also, there is a connected channel existing between the rolled-up SVs and the adjacent SSLs. The presence of the connected channel suggests that TKE can be transferred from the SSLs to the rolled-up SVs, which supports the feeding mechanism hypothesised earlier. Moving downstream, the TKE keeps increasing particularly through the canopies of the LSBs and the turbulent flow behind the LSBs remains attached to the wall (
$x=-0.32$
) until it finally separates (
$x=0$
).
Figure 18. Contour plots of time-averaged TKE (turbulent kinetic energy) for the WLE-4 case at various streamwise locations.
3.4 Influences on the aerodynamic forces
Based on the investigations made in this section, some comments are provided here in terms of the mean aerodynamic performance of the undulated aerofoils. As far as the current flow condition (deep stall) is concerned, the undulated case (WLE-8) produces a higher
$C_{L}$
and a lower
$C_{D}$
with their ratio (
$C_{L}/C_{D}$
) increased in comparison to the baseline case (SLE-8), as was the same trend in Hansen et al. (Reference Hansen, Kelso and Dally2011) – see table 2. This may be a combined outcome of the following three major events taking place in the WLE cases: (i) the appearance of a large low-pressure zone near the leading edge created by the LSBs (see figure 5); (ii) the reattachment of flow behind the LSBs resulting in a decreased size of the wake region; and, (iii) the weakening of shed vortices (to be detailed in § 4). Figure 19 shows that the size of the wake region is substantially reduced in the WLE case due to the flow reattachment taking place behind the LSB group. One of its aerodynamics consequences revealed in figure 20(a) is that the sectional lift (per unit span) estimated by calculating the circulation around the aerofoil (the Kutta–Joukowski theorem) is significantly higher in the WLE case across the majority of the span where the wake is relatively deflated. The circulation was calculated on a circle of
$R=1$
with its centre located at the centre of the aerofoil,
$(x_{0},y_{0})=(0,0)$
. The net increase in lift coefficient estimated based on the Kutta–Joukowski theorem is
$\unicode[STIX]{x1D6E5}\langle C_{L}\rangle =\langle C_{L}\rangle _{WLE}-\langle C_{L}\rangle _{SLE}=0.064$
, which precisely matches the value provided in table 2.
Figure 19. Iso-surfaces of
$\langle u\rangle =0$
in a perspective view for the (a) SLE-8 and (b) WLE-8 cases showing the size of the wake regions.
Figure 20. Relative changes in aerodynamic forces between the SLE-8 and WLE-8 cases in terms of: (a) the variation of sectional lift along the span estimated from circulation (
$\unicode[STIX]{x1D6E4}_{z}$
) around the aerofoil; and, (b,c) the piecewise distributions of lift and drag along the chord calculated from wall pressure. The sectional lift coefficient is obtained via Kutta–Joukowski theorem (
$c_{l}=2\unicode[STIX]{x1D6E4}_{z}U_{\infty }^{-1}L_{c}^{-1}$
). The relative differences between the SLE and WLE cases are denoted by
$\unicode[STIX]{x1D6E5}\{\cdot \}=\{\cdot \}_{WLE}-\{\cdot \}_{SLE}$
. Contributions from viscous stresses are excluded.
Figure 20(b,c) shows piecewise distributions of the lift and drag increase/decrease along the chord (
$\text{d}(\unicode[STIX]{x1D6E5}\langle C_{L}\rangle )/\text{d}x$
and
$\text{d}(\unicode[STIX]{x1D6E5}\langle C_{D}\rangle )/\text{d}x$
) between the SLE-8 and WLE-8 cases. It is evident in the figures that the low-pressure zone located in the LSB group is the major cause of the lift growth in the WLE case. The lift growth begins to disappear at the end of the central LSB where the flow impinges on the wall and stagnates. In the meantime, it is also shown that the low-pressure zone contributes towards increasing the drag as well. However, the overall drag reduction is achieved in the rear part of the aerofoil. It is suggested that the lower drag in the rear part is associated with the weakening of shed vortices, i.e. a loss of periodic suction pressure at the rear of the aerofoil. The loss of suction pressure seems to be contributing towards reducing the lift as well but it is relatively less pronounced, presumably because it affects the upper and lower surfaces both. The relative changes in lift and drag calculated for two split zones (front and rear) of the aerofoil are summarised in table 4.
Table 4. Differences in lift and drag forces between the WLE-8 and SLE-8 cases,
$\unicode[STIX]{x1D6E5}\langle C_{L}\rangle =\langle C_{L}\rangle _{WLE}-\langle C_{L}\rangle _{SLE}$
and
$\unicode[STIX]{x1D6E5}\langle C_{D}\rangle =\langle C_{D}\rangle _{WLE}-\langle C_{D}\rangle _{SLE}$
, over two different sections of the aerofoil chord (front and rear) where the front section is from the leading edge to the end of the central LSB (
$x\in [-0.5,-0.2]$
) and the rear section is the rest
$x\in [-0.2,0.5]$
. Contributions from viscous stresses are excluded.
4 Unsteady characteristics
In this section, the focus is moved on to the unsteady flow features of the undulated aerofoils. The effect of WLEs on the aerodynamic force fluctuations is investigated first. Secondly, the enhanced flow variety in the spanwise direction influencing the large-scale flow behaviour is discussed. Lastly, additional discussions are provided with regard to the shear layers generated at the leading-edge area.
4.1 Force fluctuations and von-Kármán vortex shedding
It is shown in table 2 that the level of lift and drag fluctuations in time (
$\unicode[STIX]{x1D70E}_{C_{L}}$
and
$\unicode[STIX]{x1D70E}_{C_{D}}$
) is reduced in the WLE cases compared to the SLE counterpart. For a stalled aerofoil, the force fluctuations are usually associated with von-Kármán vortex shedding in the rear wake. The current results also exhibit such a periodic vortex shedding at a particular frequency. In this paper, a non-dimensional frequency is defined as follows in order to identify the Strouhal number of the vortex shedding:
$$\begin{eqnarray}f^{\ast }=\frac{f\,L_{c}\sin \unicode[STIX]{x1D6FC}}{u_{\infty }}.\end{eqnarray}$$
For bluff bodies at the current Reynolds number, the Strouhal number is typically in the range of
$St\approx 0.19{-}0.20$
(Lienhard Reference Lienhard1966). In order to assess the nature of the force fluctuations of all simulations, the power spectral density (PSD) of a time signal
$q(t)$
is defined as
$S_{qq}(f)=|\hat{q}(f)|^{2}$
where
$\hat{q}(f)$
is the Fourier transform of the signal. Figure 21 shows the PSD of lift and drag fluctuations, where the SLE case exhibits a clear peak at
$f^{\ast }=0.2$
indicating von-Kármán vortex shedding as mentioned above. A secondary peak is also identified at the first harmonic of the fundamental frequency. It is envisaged that the second peak is related with the repetitive low-pressure contribution from each shed vortex, whereas the fundamental peak is related with the change of circulation by the induced velocity between two successive counter-rotating vortices.
Figure 21. PSD of
$C_{L}$
and
$C_{D}$
fluctuations for the WLE cases compared with the SLE-8 case.
In the meantime, the WLE cases (except WLE-2) show no obvious signature of the vortex shedding at
$f^{\ast }=0.2$
. The WLE-to-SLE ratios of the PSD magnitude at this frequency in figure 21 are shown in table 5 where
$\tilde{S}_{C_{L}}(f^{\ast })$
and
$\tilde{S}_{C_{D}}(f^{\ast })$
are defined as
$$\begin{eqnarray}\tilde{S}_{C_{L}}(f^{\ast })=\frac{S_{C_{L}C_{L}}(f^{\ast })|_{WLE}}{S_{C_{L}C_{L}}(f^{\ast })|_{SLE-8}}\quad \text{and}\quad \tilde{S}_{C_{D}}(f^{\ast })=\frac{S_{C_{D}C_{D}}(f^{\ast })|_{WLE}}{S_{C_{D}C_{D}}(f^{\ast })|_{SLE-8}}.\end{eqnarray}$$
The weakening of fluctuations at the shedding frequency is evident in lift and drag both. Also, the WLE cases generate weaker fluctuations in a wide range of frequencies rather than locally at the shedding frequency. It seems that the WLE cases yield vortex shedding at slightly higher frequencies close to
$f^{\ast }\sim 0.3$
. For example, in figure 22, the WLE-8 case displays more number of trailing vortices (counting those with positive
$\unicode[STIX]{x1D714}_{z}$
) than the SLE-8 case within a given distance from the trailing edge. However, the vortex shedding has almost insignificant impact on the aerodynamic force fluctuations in the WLE cases as shown in figure 21.
Figure 22. Instantaneous contour plots of (a,b) spanwise vorticity and (c,d) perturbed pressure (
$p^{\prime }=p-\langle p\rangle$
) in the rear wake region comparing the SLE-8 and WLE-8 cases. The plots are taken from T4 cross-section (denoted in figure 7).
Table 5. Relative changes in the PSD of lift and drag fluctuations for the WLE cases compared to the SLE-8 case, obtained from figure 21 at the vortex shedding frequency (
$f^{\ast }=0.2$
).
Figure 23 shows the PSD of the pressure drag (
$C_{Dp}$
) and skin-friction drag (
$C_{Df}$
) fluctuations for the SLE-8 and WLE-8 cases. Due to the fact that
$\langle C_{Dp}\rangle$
is dominant over
$\langle C_{Df}\rangle$
, their fluctuation levels also show a large difference. It is found that the WLE case yields significantly lower fluctuations of
$C_{Dp}$
at mid-to-high frequencies (figure 23
a) and of
$C_{Df}$
at low-to-mid frequencies (figure 23
b). The reduced fluctuation level around
$f^{\ast }\sim 0.2$
is consistent for both components. However, the fluctuations in
$C_{Df}$
seem to increase at high frequencies on the contrary. The small-scale near-wall turbulence in the large zone of attached flow may give rise to the high-frequency components of
$C_{Df}$
fluctuations. Nevertheless, the overall level of fluctuations is reduced in both
$C_{Dp}$
and
$C_{Df}$
for the WLE case where the reduction is more pronounced in
$C_{Df}$
as shown in table 3. It should be noted, however, that the time-averaged skin-friction drag (
$\langle C_{Df}\rangle$
) is increased in the WLE case as also shown in table 3.
Figure 23. The PSD of
$C_{Dp}$
and
$C_{Df}$
fluctuations: a comparison between the SLE-8 and WLE-8 cases.
$C_{Dp}$
and
$C_{Df}$
denote pressure and skin-friction drag coefficients, respectively.
4.2 Spanwise coherence
Figure 22 shows an evidence of the vortex shedding becoming significantly weaker in the WLE cases. It is shown that the spanwise vortical structures downstream of the SLE case remain strong even after travelling several chord lengths. Those in the WLE case are much weaker albeit still existing. It is also shown in the figure that the SLE case produces a large coherent vortical structure rolling up (anticlockwise) at the trailing edge, which does not appear in the WLE case. The deteriorated vortex shedding in the WLE case can be explained by the loss of coherence of the flow in the spanwise direction associated with the leading-edge modifications. A magnitude-squared coherence spectrum
$C_{p_{1}p_{2}}$
of pressure fluctuations between two different locations
$p_{1}(t)$
and
$p_{2}(t)$
is defined as
$$\begin{eqnarray}C_{p_{1}p_{2}}(f)=\frac{|S_{p_{1}p_{2}}(f)|^{2}}{S_{p_{1}p_{1}}(f)S_{p_{2}p_{2}}(f)},\end{eqnarray}$$
where the two-point cross-spectrum and the cross-correlation of the signals are given by
$$\begin{eqnarray}S_{p_{1}p_{2}}(f)=\int _{-\infty }^{\infty }R_{p_{1}p_{2}}(t)\text{e}^{-2\unicode[STIX]{x03C0}\text{i}ft}\,\text{d}t\quad \text{and}\quad R_{p_{1}p_{2}}(t)=\langle p_{1}(\unicode[STIX]{x1D70F})p_{2}(\unicode[STIX]{x1D70F}+t)\rangle .\end{eqnarray}$$
The spanwise coherence is computed between two points along a straight line in the
$z$
-direction with one of them (
$p_{1}$
) fixed as a reference. This process is repeated for 31 different streamwise locations along the chord. Figure 24 reveals the signature of the spanwise coherence at the shedding frequency
$f^{\ast }=0.2$
for the SLE-4 and WLE-4 cases obtained on a horizontal plane of
$y=0.15$
. It is clear in the figure that the WLE case shows much lower spanwise coherence signature at the vortex shedding frequency, particularly in the region where the group of LSBs is located. This suggests that the highly promoted three-dimensionality of the flow near the leading edge discourages the formation of strong coherent vortical structures downstream.
Figure 24. Contour plots of the magnitude-squared coherence of pressure fluctuations between two points along the span obtained at the frequency of
$f^{\ast }=0.2$
on a horizontal plane of
$y=0.15$
, for the (a) SLE-4 and (b) WLE-4 cases. The location of the reference signal (
$p_{1}$
) is denoted by the thick dashed line.
4.3 Leading-edge vortex dynamics
Returning our focus back to the leading-edge area, the unsteady characteristics of the long spanwise coherent structures shown earlier in figure 12 are investigated here. The long spanwise coherent structures are formed in the leading-edge shear layers due to a Kelvin–Helmholtz instability (Watmuff Reference Watmuff1999; Yarusevych, Sullivan & Kawall Reference Yarusevych, Sullivan and Kawall2009). It has been discussed in § 3.2 that they tend to turn towards the streamwise direction in the WLE cases and interact with the LSBs nearby.
A series of pressure probes are distributed along a curve where the root-mean-square of the pressure fluctuations is highest (tracing the shear layer), on two different
$xy$
-planes for comparison. The PSDs of the pressure fluctuations obtained are shown in figures 25 and 26 for the SLE-2 and WLE-4 cases, respectively. In the SLE-2 case (figure 25), there is no major difference between the cross-sections T1 and T2 and they display strong local peaks at
$f^{\ast }\approx 6{-}6.5$
across the first two probe points indicative of periodic vortex shedding. As the probe moves downstream, the signature of the periodic mode fades away and the overall energy level rises until it reaches a saturation. In the WLE-4 case, however, there are significant differences between the cross-sections T2 and T4 due to the enhanced flow three-dimensionality. In figure 26(b) obtained from the cross-section T2 cutting through the central LSB, a local peak is observed at
$f^{\ast }\approx 7{-}7.5$
at the two most upstream probe points. The peak broadens out as the probe moves downstream and the energy seems to be transferred to lower frequencies. The overall energy level also decreases with the probe moving downstream. In contrast, figure 26(c) obtained from the cross-section T4 cutting through the SSL shows a major peak at a much lower frequency
$f^{\ast }\approx 3$
in the upstream region. In this case, the overall energy level does not seem to decay as much as that from the cross-section T2 although the energy cascade towards lower frequencies is still apparent.
Figure 25. PSD of fluctuating pressure in the SLE-2 case obtained at some probe points: (a) the locations of the probe points aligned on the shear layer; (b) PSDs obtained at the probe points on T1 cross-section; and, (c) those on T2 cross-section. The positions of T1 and T2 are denoted in figure 7.
Figure 26. PSD of fluctuating pressure in the WLE-4 case obtained at some probe points: (a) the locations of the probe points aligned on the shear layer; (b) PSDs obtained at the probe points on T2 cross-section; and, (c) those on T4 cross-section. The positions of T2 and T4 are denoted in figure 7.
The vortex dynamics taking place on the cross-section T4 (in relation with figure 26 c) is displayed in figure 27. In figure 27(a), the first two small vortices A and B eventually merge into bigger ones (C and D). In figure 27(b), A and B have already merged into AB while C and D keep drifting without noticeable dissipation. Figure 27(c,d) shows a slightly different events on the same cross-section. In this case, two triads of vortices (A-B-C and D-E-F in figure 27 c) promptly merge together into ABC and DEF (figure 27 d). However, in this sequence, significant dissipation of the vortices is exhibited further downstream compared to the previous case. This process of vortex merging and dissipation taking place intermittently seems to result in the wide broadband spectra shown in figure 26(c).
Figure 27. Instantaneous contour plots of pressure coefficient in the WLE-4 case, obtained on the cross-section T4 intersecting the SSL region. The position of T4 is denoted in figure 7.
In contrast, the vortical structures viewed on the cross-section T2 that cuts through the central LSB seem to undergo a much more rapid merging process as shown in figure 28. In figure 28(a), the small vortices A, B and C can be individually recognised. However, a short time later (figure 28 b), B and C have already merged into BC. The continued merging of A–BC into ABC is complete in figure 28(c). In figure 28(d), the vortical structures start to dissipate in the turbulent boundary layer downstream of the LSB. Since the vortical structures are restricted to a much smaller space to roll up/down around each other due to the presence of the wall, the merging process starts almost immediately after shedding. The high shear near the wall also promotes a much rapid dissipation of the vortical structures due to straining and stretching in the streamwise direction. The combination of the fast merging and dissipation process could explain the widening and disappearance of the spectral peak observed in figure 26(b).
Figure 28. Instantaneous contour plots of pressure coefficient in the WLE-4 case, obtained on the cross-section T2 intersecting the central LSB. The position of T2 is denoted in figure 7.
5 Concluding remarks
The flow around a deep-stalled NACA0021 aerofoil with a wavy leading edge (WLE) has been investigated by using high-order accurate large-eddy simulations. The current study presented comprehensive analyses on both time-averaged and unsteady results. In the WLE cases, low-pressure spots behind trough areas were identified, which were consistent with the previous observations in the literature. The simulation data showed that the low-pressure spots were in fact the signature of LSBs. It was revealed in this study that these LSBs were not present after every trough of the WLEs and the distribution of the LSBs varied with the number of wavelengths (spanwise domain size) employed in the simulations. More importantly, there was a consistent trend that the LSBs formed a group together and a large zone of attached flow was developed behind the LSB group. The size of the LSB group and the attached flow zone grew as the number of wavelengths used was increased from two to eight. As far as the current simulations are concerned, there was always one trough area with an SSL occurring almost immediately after the leading edge regardless of the number of wavelengths used.
It was found that, in the WLE cases, a pair of prominent streamwise vortices were created due to two opposite vorticity layers overlapping and rolling up against each other. It was suggested that the rolled-up streamwise vortices acted as a buffer to prevent the adjacent SSLs from penetrating into the LSB group. It was also suggested that the longevity of the rolled-up streamwise vortices was maintained by the feeding of streamwise vorticity and turbulent kinetic energy from the adjacent SSLs through a connected channel. Due to the feeding mechanism, only two outermost vortices in contact with the SSLs were sustained and lasted until the SSLs were dissipated away. This allowed for the appearance of the attached flow zone behind the LSB group. The variety of flow features (three-dimensionality) in the WLE cases resulted in a reduced spanwise coherence of the flow, and contributed to producing a significantly weaker vortex shedding at the rear of the aerofoil compared to the SLE counterpart.
The current result agreed well with the previous experimental measurement by Hansen et al. (Reference Hansen, Kelso and Dally2011) in terms of time-averaged aerodynamic performance. The growth in lift and reduction in drag with the WLE geometry were confirmed. Some explanations to the enhanced aerodynamic performance were provided in this paper. The reduced size of wake due to the attached flow region resulted in a higher circulation and therefore lift. Locally, the low-pressure zone created under the LSB group was the major contribution to the growth in lift, although it also contributed towards increasing the drag. The reduction in drag was achieved mainly by the weakening of shed vortices (loss of suction) at the rear of the aerofoil. It is worth noting that the skin-friction drag was higher in the WLE cases due to the attached flow behind the LSB group although this was outweighed by the pressure drag. The unsteady aerodynamics and flow characteristics have also been investigated. The reduced strength of the shed vortices in the WLE cases resulted in significantly lower levels of lift and drag fluctuations than those of the SLE case, at both the peak and broadband frequencies alike. In addition, the rapid broadening and decaying of the turbulent energy spectra in the leading-edge vortex dynamics were evident with the WLE geometry.
Acknowledgements
The authors gratefully acknowledge high-performance computing facilities and services of the UK National Supercomputer ARCHER via the support of UK Turbulence Consortium (EP/L000261/1) and the local IRIDIS4 at the University of Southampton in the completion of this work. All data supporting this study are openly available from the University of Southampton repository at https://doi.org/10.5258/SOTON/404212.
