Direct numerical simulations of the flow around wings with spanwise waviness at a very low Reynolds number

Inspired by the pectoral flippers of the humpback whale, the use of spanwise waviness in the leading edge has been considered in the literature as a possible way of improving the aerodynamic performance of wings. In this paper, we present an investigation based on direct numerical simulations of the flow around infinite wavy wings with a NACA0012 profile, at a Reynolds number Re = 10 0 0 . The simulations were carried out using the Spectral/hp Element Method, with a coordinate system transformation employed to treat the waviness of the wing. Several combinations of wavelength and amplitude were considered, showing that for this value of Re the waviness leads to a reduction in the lift-to-drag ratio ( L / D ), associated with a suppression of the fluctuating lift coefficient. These changes are associated with a regime where the flow remains attached behind the peaks of the leading edge while there are distinct regions of flow separation behind the troughs, and a physical mechanism explaining this behaviour is


Introduction
The possibility of using wings with a wavy leading edge as a way of obtaining improved aerodynamic performance started receiving attention after the work of Fish and Battle [3] , where the morphology of the pectoral flipper of the humpback whale ( Megaptera novaeangliae ) was analysed with a focus on its hydrodynamic performance. These flippers have protuberances on the leading edge, and this was suggested to act as a mechanism to delay the stall, allowing the flipper to maintain a high lift coefficient at high angles of attack, giving the whale a good maneuverability.
The idea that leading edge protuberances could delay stall gained support with the work of Miklosovic et al. [9,10] . They presented experiments for full-span and half-span wings with a NACA0020 profile in configurations with and without leading edge waviness. For the half-span model, which had a planform similar to the flipper of the humpback whale, the Reynolds number was around 6 × 10 5 and the modified wing led to an increase in the stall angle. This increase in the stall angle contributed to an increase in the maximum lift coefficient of the wing. However, for the full-span model, for which the Reynolds number was around 2.7 × 10 5 , their results show that the protuberances lead to a premature stall, being beneficial only in the post-stall regime. Also, the experiments for full-span wings presented in [5] showed the same behaviour, with the modified leading edge causing a premature stall.
Although these first studies about the effect of leading edge protuberances showed a strong distinction between the behaviour of full-span and half-span models, more recent works suggest that the main factor affecting the results is the Reynolds number. First, Stanway [14] presented experiments for a model similar to the half-span wing of Miklosovic et al. [10] , but considering different values of Re between 4 × 10 4 and 1.2 × 10 5 . Only for the highest value of Re considered the waviness caused an increase in the maximum lift coefficient, indicating that the value of Re has an important role in determining whether the use of wavy leading edges will improve aerodynamic performance. Another study which supports the importance of the Reynolds number effect on this flow is that of Hansen et al. [4] . They performed experiments with rectangular wing mounted in both full-span and half-span configurations, in an attempt to isolate the influence of the wing tip on the results. The effect of using a wavy leading edge was similar in both cases, indicating that three-dimensional effects related to the wing-tip have a secondary importance in the effectiveness of the waviness.
Despite the significant number of studies on this problem, a definite explanation to how the leading edge waviness affects the flow is still absent. In this paper, an extensive study based on direct numerical simulations is presented for the flow around infinite wavy wings with a NACA0012 profile at Re = 10 0 0 . Although this Reynolds number is much lower than most practical applications, it is our belief that the conclusions presented here can help in the understanding of the mechanisms responsible for the behaviour observed at higher Reynolds numbers.
The paper is organized as follows. Section 2 briefly presents the numerical methods and describes the setup employed in the simulations. Then, Section 3 presents the results, and finally Section 4 contains the conclusions of this work.

Numerical methods
We consider an incompressible viscous flow, which is governed by the Navier-Stokes equations. Assuming a unit density, these equations can be written as: where u is the velocity, p is the pressure, and ν is the kinematic viscosity. Taking the chord c as reference length and the freestream velocity U ∞ as the reference velocity, the Reynolds number is defined as Re = cU ∞ ν . The waviness of the wing was treated by a coordinate system transformation using the formulation proposed in [12] , so that the wavy geometries were mapped into the straight wing. Then, the equations were discretized following the Spectral/hp method presented in [8] , with the span direction discretized using a Fourier expansion, as proposed by Karniadakis [6] . The use of a Fourier expansion in the third direction is an efficient way of studying an infinite wing with periodic boundary conditions, and was only possible because of the simplification in the geometry provided by the coordinate transformation. The equations were then solved by time integration using the stiffly stable splitting scheme described by Karniadakis et al. [7] .

Simulations setup
All numerical simulations presented here were performed for Re = 10 0 0 , using a modified version of the incompressible Navier-Stokes solver encapsulated within the spectral/hp element code Nektar ++ [1] . The geometries are based on a NACA 0012 profile, with a small modification to obtain a zero-thickness trailing edge, since the small thickness of the trailing edge in the original profile would lead to meshing challenges. The wavy geometries were obtained by applying the following coordinate transformation to the straight infinite wing: where h is the waviness peak-to-peak amplitude, λ is its wavelength, x is the physical coordinate in the chord direction, and x and z are the chord wise and span wise coordinates in the computational domain. Fig. 1 shows an example of a geometry obtained through this transformation, identifying the waviness parameters and the orientation and origin of the coordinate system. Note that this transformation deforms both the leading and the trailing edges, and has no effect on the chord. This type of transformation was preferred because it leads to a simpler procedure when solving the equations with the mapping [12] . However, in Section 3.5 we relax this restriction by considering the effect of  waviness in the leading edge and in the trailing edge independently. As will be seen, the transformation from Eq. (2) leads to results that are equivalent to deforming only the leading edge, as is usually the case in the literature. We considered the case without any waviness (referred to as baseline) and wavy geometries with nine different combinations of the parameters λ and h . These cases are summarized in Table 1 , where a naming convention is also introduced. The parameters were chosen taking into account the range of parameters encountered in the literature [2,4] , with the amplitudes adjusted according to the wavelength in order to have the same ratios h λ for each wavelength.For each case, simulations were performed for angles of attack between α = 0 • and α = 21 • , at increments of α = 3 • .
For all simulations the spatial discretization in the xy plane consisted of a mesh with 550 quadrilateral elements, with the solution represented by 8th order polynomials inside each element. This mesh, which is shown in Fig. 2 , extends from −10 c to 10 c in the chord direction x , and from −15 c to 15 c in y . The z di- rection was represented by a Fourier expansion with 16 degrees of freedom extending a span of one chord. Each simulation was computed for 100 time units with a time-step of 5 × 10 −4 , and the last 10 time units were considered to obtain the time-averaged forces. When lift coefficient fluctuations are reported, these were obtained by extending the simulation for an additional 100 time units, since these results were more sensitive to the averaging interval, especially at high angles of attack. These parameters were chosen based on an extensive convergence study, which consisted of systematically varying the relevant parameters (polynomial order, mesh size, time step and number of Fourier modes) independently, and choosing a value which provides a sufficiently small difference com pared with the most refined configuration. For tests with α = 5 • , increasing the polynomial order from 8 to 14 or doubling the number of Fourier modes we observe a change on the mean forces lower than 0.2%. This same threshold is respected if we consider the effect of using a lower time-step of 2 . 5 × 10 −4 , or using a larger domain extending from −15 c to 10 c in x and from −20 c to 20 c in y .
In all simulations, the boundary conditions for the velocity consisted of: unit velocity oriented to provided the appropriate angle of attack on the left, bottom and top boundaries; zero velocity on the wing surface; and zero normal direction velocity derivative on the right boundary. The pressure was set to zero in the outflow, while in the inflow and wall boundaries a high-order pressure boundary condition was employed [7,12] . Periodic boundary conditions are automatically imposed along the span, due to the use of a Fourier expansion in the z direction.  Fig. 5 . For the baseline wing, the forces have an oscillatory behaviour which can be associated to vortex shedding at this high angle of attack. The L05h05 geometry shoes a similar behaviour, but with smaller fluctuations, while in the larger amplitude case L05h10 these fluctuations are suppressed, with C L corresponding only to small fluctuations with a lower frequency.

Separation and recirculation zones
In order to obtain a better understanding of the characteristics of the flow which lead to the behaviour described above, Fig. 6 shows skin friction lines on the wall of the wing for the baseline, L05h05 and L05h10 geometries. The skin friction lines represent the orientation of the projection on the wing surface of the shear force (the shear stress tensor dotted with the surface normal vector). This figure represents top views for α = 9 • , α = 12 • and α = 15 • , with the flow oriented from left to right. Also, to make visualization easier, the colours represent the orientation of the skin friction in the chord direction, with blue regions corresponding to reversed flow associated with boundary layer separation. For the baseline wing, the flow is approximately two-dimensional in the range of angles of attack considered in the figure, with an increase in the angle of attack leading to flow separation closer to the leading edge. In the other cases, the use of the waviness tends to prevent separation in the regions corresponding to the peaks of the To further illustrate the behaviour discussed above using skin friction lines, Fig. 7 shows the recirculation regions of one the cases of Fig. 6 , namely the case L05h10 with α = 12 • . This visualization was obtained considering the regions where the stream- wise velocity is negative. From this figure, it is clear that for this case, the boundary layer separation is restricted to the waviness troughs.

Surface pressure distribution
A good starting point for understanding the previous results is looking at the pressure distribution on the wing surface, as it is related to flow separation and to the aerodynamic forces. Fig. 8 (a) shows the pressure coefficient distribution on the wing surface for the case L05h10 with α = 12 • along two different cross-sections, one at the waviness peak and the other on the valley, and also the distribution on the baseline wing as a reference. From this plot, it is clear that there is a strong reduction in the leading edge suction for the section on the waviness peak. The section containing the waviness valley also experiences a reduction in the suction peak, although it is much less severe. The reduction in the aerodynamic forces observed earlier can be attributed to this loss of suction on the upper surface of wing, as becomes clear when we look at the areas inside the −C p curves. Thus, even though there is a reduction in the separation region, the loss of leading edge suction prevents this from being converted to an increase of lift.
Another important result from the pressure distribution of Fig. 8 (a) is that, as the flow progresses along the chord, there is a tendency for the pressure to equalize along the sections, eliminating spanwise pressure gradients. These gradients are almost entirely eliminated at x = 0 . 4 , and since the suction peak on the valley section occur behind the corresponding point on the peak section, the former has a shorter distance to recover the pressure when compared to the latter. This, combined with the stronger suction in the valley makes this section experience stronger adverse pressure gradients, explaining why separation is restricted to these parts of the wing. Fig. 8 (b) presents the (streamwise) tangential pressure gradients corresponding to the previous pressure distribution. While the gradients on the valley section are comparable to the baseline case, the peak section faces weaker gradients, confirming the previous argument.
Also, it should be emphasized that, at least for the low Reynolds number considered here, the fact that the valley section has a shorter distance to recover the pressure is related to the elimination of spanwise pressure gradients, and not to a physical restriction imposed by a shorter chord as has been suggested in the literature, for example in [15] . In fact, in the simulations corresponding to Fig. 8 the chord is constant along the span, and therefore there is no chord length effect present.

Physical mechanism
The previous sections presented results showing that, for moderate angles of attack, the presence of waviness in the wing leads to a reduction of L / D , which is associated with a flow regime where separation is restricted to the regions behind the waviness valleys. Also, this can be explained by considering the pressure distribution around the wing surface. Therefore, the only remaining point to obtain a good understanding of the flow is explaining how this pressure distribution is formed. This section addresses this issue, proposing a physical description which can justify the observed behaviour, based on a simple sectional argument for the region close to the leading edge.
The main point to understanding how the waviness affects the flow is to note that, by deforming the wing, spanwise pressure gradients are created. To illustrate this, Fig. 9 (a) shows the spanwise pressure gradients contours at z = 0 . 125 , which is a plane halfway between a peak and a valley of the waviness, for the case L05h10 with α = 12 • . The closest peak is at z = 0 . 0 and the closest valley at z = 0 . 25 , and thus a positive gradient is oriented from peak to valley and a negative gradient from valley to peak. Close to the leading edge, where x < 0, the flow has already reached the wing in the peak of the waviness, and therefore the pressure gradient is negative on the lower surface of the wing (due to the stagnation point) and positive on the upper surface of the wing (due to the suction peak). This pressure gradient accelerates the flow in the span direction, generating a spanwise flow moving away from the protuberance on the lower surface and towards it on the upper surface, as shown in Fig. 9 (b). This behaviour can be seen as an attempt of the flow to move around the protuberance, as illustrated in Fig. 10 . In this figure, the coloured contours represent spanwise velocity at the plane x = −0 . 03 , which is close to the location of the leading edge in the peak of the waviness at x = −0 . 05 , showing how the flow moves away from the waviness peak in the lower portion of the wing, while it moves towards it in the upper part. Also, the grayscale represents the pressure on the wing surface, with the darker areas corresponding to a lower pressure (a stronger suction). The effect of the flow towards the suction peak is to increase the pressure at that point, causing the reduction of −C p at this section as highlighted in Fig. 8 . Also, this spanwise flow around the leading edge protuberance is somewhat compatible with the behaviour of the streamlines obtained by Skillen et al. [13] at a higher value of Re . However, although they observed the flow being deflected by the lower surface of the leading edge, their results do not show any flow towards the suction peak in the upper surface. They suggest that this deflection lead to the acceleration of the flow behind the troughs, generating an improved suction peak. However, this effect is not observed in our results.
The previous argument can also be interpreted in a more intuitive manner. By deforming the leading edge along the span with displacements in the chordwise direction, we hinder the wing's ability to force the flow to accelerate around the leading edge, specially in the peak section. Although this is not optimal in the sense of generating the maximum possible suction peak (and consequently the maximum possible lift), it prevents the flow from separating.
If we consider now the flow at downstream positions, we note that as the flow reaches the rest of the wing, at x > 0, the spanwise pressure gradients are reversed. The spanwise flow then acts to eliminate the spanwise pressure gradients, leading to the behaviour for the surface pressure described in the previous session. Clearly, the existence of gradients of the w velocity is only possible in the presence of streamwise vortices. However, from the contours of streamwise vorticity of Fig. 9 (c) we notice that these vortices are lifted away from the wing. Therefore, they are not expected to have a significant impact on the observed aerodynamic behaviour, being rather simply a consequence of the spanwise flow.
In conclusion, this section has used simple arguments based on the spanwise flow to explain how the pressure distribution is affected by the use of waviness. The modification of the pressure distribution leads to weaker suction peaks, necessarily leading to a more benign pressure gradient behind the waviness peak. These factors are responsible for the observed separation patterns and for the loss of lift caused by the waviness.

Effect of waviness with different shapes
As noted previously, the simulations presented so far considered a wing with uniform deformation in each xy plane, so that both the leading and trailing edges are deformed. A few simula-tions were performed in order to determine whether this is significantly different from the typical geometry presented in the literature, where only the leading edge is deformed. These simulations consisted in modulating the function ξ ( z ) from Eq. (2) in the chord direction, obtaining wings with modifications to only one of the extremities. Using this transformation, the absolute thickness of the wing is constant throughout the span, and therefore the waviness troughs have a slightly thicker profile than the baseline wing (due to the shorter chord), while the waviness peaks have a slightly thinner profile (due to the longer chord). Fig. 11 presents the lift-to-drag ratio for the case L05h10 in these different configurations. The suffix LE indicates that only the leading edge is deformed, and similarly TE for the trailing edge. It is clear that deforming the whole wing is almost equivalent to deforming only the leading edge, and therefore our results from the previous sections are comparable to the literature. Also, changes to the trailing edge have little effect on the results, at least for this Re .
Another distinction that can be made for the case L05h10_LE is that, since the chord changes along the span, we have to choose between maintaining a constant absolute thickness or a constant profile (which requires a constant relative thickness). Tests were performed for these two possibilities, but the results, which are not presented here, were almost identical in both cases.

Discussion and conclusions
We have presented a study of the effect of spanwise waviness on the flow around infinite wings at a very low Reynolds number. For moderate angles of attack, the waviness leads to decreases in both drag and lift forces, which result in a decrease in the liftto-drag ratio. There is also a significant suppression of the fluctuations in the lift coefficient. This behaviour is related to a flow regime where flow separation is restricted to the regions behind the waviness valleys, with separation behind the peaks being suppressed.
A physical explanation for the results is proposed, based on the spanwise flow induced by the waviness. This explanation leads to the conclusion that by deforming the wing, we decrease its ability to force the flow to accelerate around the leading edge, resulting in a weakening of the suction peak. The flow does not separate behind the peaks because this phenomenon is stronger in this region, and therefore it does not contain strong adverse pressure gradients.
We can obtain further insight into this flow if we consider, as a rough approximation, that the effect on the lift of weakening the  suction peak is similar to the effect of reducing the angle of attack. In the present case, the lift coefficient increases monotonically with the angle of attack. Therefore, the lower effective angle of attack caused by the waviness is expected to lead to reductions in the lift, as was in fact observed in our results. However, for higher Re where a sharp stall is observed, the stall may be delayed in portions of the wing, leading to a loss of lift before stall and an increase of lift in the post-stall, as observed in the literature. This explanation is also consistent with the smoothing of the stall observed in the literature: since the effects of the waviness are not uniform along the span, a gradual progression of the stall is expected.
As a final note, we address the issue of the Reynolds number considered here being much lower than what is expected in practice, what is in part due to the limitations associated with the high computational cost of direct numerical simulations. Ongoing research at larger Reynolds numbers show that the physical mechanisms discussed in the present paper are still present even at Re = 50 , 0 0 0 [11] , demonstrating the importance of these simulations at Re = 10 0 0 to understanding the effect of the waviness on the flow.