Control of low Reynolds number flow around an airfoil using periodic surface morphing: a numerical study

The principal aim of this paper is to use Direct Numerical Simulations (DNS) to explain the mechanisms that allow periodic surface morphing to improve the aerodynamic performance of an airfoil. The work focuses on a NACA-4415 airfoil at Reynolds number Rec = 5×104 and 0◦ angle of attack. At these flow conditions, the boundary layer separates at x/c = 0.42, remains laminar until x/c ≈ 0.8, and then transitions to turbulence. Vortices are formed in the separating shear layer at a characteristic Kelvin-Helmholtz frequency of Sts = 4.9, which compares well with corresponding experiments. These are then shed into the wake. Turbulent reattachment does not occur because the region of high turbulent kinetic energy (and therefore mixing) is located too far downstream and too far away from the airfoil surface to influence the near-wall flow. The effect of three actuation frequencies is examined by performing the simulations are on a computational domain that deforms periodically. It is found that by amplifying the Kelvin-Helmholtz instability mechanism, Large Spanwise Coherent structures are forced to form and retain their coherence for a large part of the actuation cycle. Following their formation, these structures entrain high momentum fluid into the near-wall flow, leading to almost complete elimination of the recirculation zone. The instantaneous and phase averaged characteristics of these structures are analysed and the vortex coherence is related to the phase ∗Corresponding author Email address: g.papadakis@ic.ac.uk (George Papadakis) Preprint submitted to Journal of Fluids and Structures September 29, 2017 of actuation. In order to further clarify the process of reduction in the size of recirculation zone, simulations start from the fully-developed uncontrolled flow and continue for 25 actuation cycles. The results indicate that the modification of airfoil characteristics is a gradual process. As the number of cycles increases and the coherent vortices are repeatedly formed and propagate downstream, they entrain momentum, thereby modifying the near wall region. During this transient period, the separated shear layer approaches the airfoil surface and the size of recirculation region decreases. It takes at least 15 cycles for the flow to develop a repeatable, periodic pattern.


Introduction
Civilian and military interest in Micro and Unmanned Air Vehicles (MAVs and UAVs respectively) has increased significantly in recent years. These vehicles can carry very small sensors, video cameras, listening devices and control hardware systems and they are capable of complex missions ranging from 5 surveillance and communication relay to detection of biological, chemical, or nuclear materials [1]. The combination of small length scales and low velocities results in an operating flight regime of Re c in the range of 15,000 to 500,000.
At such Reynolds numbers, a laminar boundary layer forms on an airfoil's upper surface and persists beyond the suction peak and into the pressure re-10 covery region, where it encounters an adverse pressure gradient. Viscous effects close to the wall slow down a fluid element, thereby reducing its kinetic energy.
Turbulent boundary layers can entrain higher momentum fluid from the edge of the boundary layer and mix it effectively with lower momentum fluid close to the surface, thus energizing the near-wall region. In laminar boundary layers 15 this mechanism is much slower and the energy of the near-wall flow remains very low. This makes them incapable of overcoming even modest adverse pressure gradients and prone to separation even at low angles of attack, resulting in high drag and loss of lift [2, 3,4,5,6].
Flow control can counter such unfavorable behavior and potentially lead to 20 considerable performance improvements. Given an imposed pressure field, the principle strategy in separation control is to add momentum to the very nearwall region. Passive flow control techniques, like Vortex Generators (VG) [7,8,9] or surface roughness [10,3,11] are easily implemented, but their performance is optimal only at design conditions. An aircraft would encounter a variety of 25 flow conditions during a single mission, thus it is doubtful whether such control strategies can have a net positive effect in practice.
In contrast, active control approaches remove the drawback associated with passive control at off-design conditions. This form of control has been associated mainly with periodic injection or suction of fluid from a boundary layer. 30 Periodic forcing was first used by Schubauer and Skramstad [12] who introduced periodic perturbations in a laminar boundary layer to trigger a known instability (the Tolmien-Schlichting wave). More recently, the most widely researched periodic active flow control methods are Synthetic Jets (or Zero-Net-Mass-Flux jets) [13,14,15,16] and Direct-Barrier-Discharge (DBD) Plasma actuators [17]. 35 The formation of Large Coherent Structures (LCS) in the separated shear layer has been widely observed [13,18] and the periodic actuation is believed to be capable of accelerating and regulating their production. These structures are the essential 'building blocks' of the mixing layer [11] and are responsible for transporting momentum across it. The triggering of LCS is therefore an efficient 40 method for control of mixing in the separated shear layer and consequently, periodic actuation can increase the transfer of high momentum fluid, enhancing entrainment. This principle forms the basis of separation control by periodic excitation.
Since the effectiveness of the method is largely determined by the receptivity 45 of the flow to the imposed disturbances, these have to be of the right scale and be introduced at the right location. The control authority of periodic excitation has been found to exhibit a highly non-monotonic variation with actuation frequency which suggests the presence of rich flow physics [13]. Benard and Moreau [19] investigated experimentally the use of DBD plasma actuators for separation 50 control. When the actuator was activated in an initially quiescent flow field, it produced a large counter-clockwise rotating vortex. However, no such vortex was observed when used in the actual flow, but instead an amplification of a clockwise rotating LCS was found. The velocity of the gas produced by the actuator alone was one order of magnitude lower than the external velocity, but 55 still produced a dramatic effect. In agreement with what was stated above, this led the authors to conclude that the momentum transfer did not come from the actuator directly but the actuator instead acted as a catalyser that reinforced an already existing instability.
Many numerical simulations have been reported in the literature to elucidate 60 the effect of periodic actuation on the separation flow characteristics. Postl et al. [20] used Direct Numerical Simulations to investigate pulsed vortex generator jets (pVGJ) for separation control. They also observed the production of LCS and were able to confirm their spanwise coherence. The control mechanism was related to the formation of these spanwise LCS. The authors argued that pVGJ are able to generate spanwise LCS as a result of the hydrodynamic instability mechanisms of the separated shear layer. Since an inflection point is present in the velocity profile, the flow is susceptible to the Kelvin-Helmholtz (K-H) instability. More recent studies by Marxen et al. [21], Sato et al. [17] and Buchmann et al. [14] also observe spanwise LCS that are created due to the K-H 70 instability mechanism. When steady forcing and pulsed forcing were compared, Postl et al. [20] found that pulsed blowing is significantly more effective when the same momentum coefficient was used for the actuation.
Despite the improved efficiency of periodic-based active control methods with respect to both steady active and passive control, they have been diffi-75 cult to apply in practice owing to the complexity of the control devices [17].
The review paper of Cattafesta and Sheplak [22] provides many details on the advantages and disadvantages of various actuators. The development of new actuation devices and material systems has enabled novel approaches to periodic flow control to be explored. Munday et al. [23] used a thin, flexible piezoelectric 80 THUNDER actuator, developed at NASA, to morph the surface of an airfoil.
When embedded in a surface or attached to flexible structures such actuators provide a distributed force with little power consumption. They are also very light and easy to integrate to the surface of an airfoil, thus enabling their practical use and maximizing the possible aerodynamic gains. Munday et al. [23] 85 performed both static and dynamic morphing tests. While their static tests did not prove very successful, dynamic actuation was found to significantly reduce flow separation. Local wall oscillation in the suction side very close to the leading edge was also considered in the 2D numerical study of Kang et al. [24]. The authors found that when the frequency of oscillation locks into the frequency 90 of the flow, vortices are formed close to the leading edge, they roll along the surface and maintain a low pressure distribution, thereby enhancing lift. The Reynolds number was very low, equal to 5000, and the dominant flow frequency was due to vortex shedding.
The effectiveness of periodic surface morphing of the suction side has been 95 recently demonstrated experimentally for a NACA-4415 airfoil at Re c = 50, 000 and various angles of attack [25,26]. In this paper, we use DNS to investigate this form of actuation for the same airfoil and Reynolds number. Due to the cost of computations, we restrict our simulations to 0 • angle of attack. The surface morphing poses additional computational challenges because it necessitates the 100 use of deforming meshes. The central aim is to elucidate further the effect of this form of actuation on the flow characteristics.
The paper is organised as follows. Section 2 presents the computational methodology and discretistion details. This is followed by section 3 with results from the static (i.e. unactuated) airfoil, while section 4 presents results from 105 3 actuation frequencies. Finally section 5 summarises the main findings of the paper.

Governing equations and solution details
Due to the variation of the solution domain with time when the airfoil suction side is periodically morphed, the integral form of the incompressible Navier-Stokes and continuity equations is considered (refer to [27] for details about their derivation). For an arbitrarily moving and deforming control volume (CV) with bounding surface (S), the equations take the form where u i , u pi (i = 1, 2, 3) are the velocity components of the fluid and the control 110 volume surface respectively in the Cartesian directions  The finite volume method was used to discretise these equations. In order to avoid the generation of artificial mass sources (or sinks) when the cell faces move relative to each other, the discretisation should satisfy the Space Conservation Law (refer to [28,27]) This law can be thought of as the continuity equation at vanishing fluid velocity.
The set of equations (1) and (2) was discretised and solved using the OPEN-  In the DNS study by Jones et al. [30], for flow past a NACA-0012 airfoil at 150 Re c = 5×10 4 and 5 • angle of attack, the boundaries were located at a distance 5.3c. This modest domain size was found to be sufficient to capture the outer potential flow (computations in a larger domain with radius 7.3C showed small variation in the pressure distribution). In the present simulations, the angle of attack is 0 • , so the lateral extent of pressure disturbance due to the presence of 155 the airfoil is limited. A domain size of 10c is considered therefore sufficient.
The grid boundary node distribution is specified on the airfoil surface and the outer boundary. In total 1024 points were specified in the ζ direction around the airfoil surface, 220 in the η direction and 1084 in the wake. The domain was extruded 0.2c in the z direction to form a 3D grid. This has been found to 160 be sufficient for capturing the 3D flow behaviour in previous studies [30,31]. A coarse grid with 16 planes in the z direction is used to initialize the simulations.
A finer grid, with 64 planes in the z direction, is then employed for better resolution. The total number of cells for the fine grid is 45 million. The internal mesh is generated using linear transfinite interpolation between 165 the corresponding boundary coordinates on the airfoil and the outer boundary [32]. An elliptic refinement with a Neumann boundary condition is implemented to increase grid orthogonality at the airfoil surface and improve accuracy in the near-wall region. The resulting near-wall grid lines in the vicinity of the leading and trailing edges can be seen in Figure 2. Using the elliptic-Neumann 170 refinement technique the maximum non-orthogonality is reduced from 73 • to 26 • , and the mean from 19 • to 6 • .
A pressure outlet was used at the outflow exit (shown as red line in Figure 1) and a uniform free stream velocity was used on the rest of the external boundary (shown as blue curve in Figure 1). The no-slip condition was enforced on the were placed in the same region).

Implementation of periodic surface morphing 200
In the experimental investigation of surface morphing [25,26], the flow was found to have little effect on the surface motion, thereby in the simulations it was not deemed necessary to solve a coupled fluid-structure interaction problem.
The amplitude distribution and the angular frequency define completely the motion of the moving boundary. In the physical model, the leading and trailing 205 edges of the deforming portion of the surface -herein referred to as the control surface -are attached to the airfoil body by a rigid bond and a slide joint respectively. As a result, the wall-normal displacement of the control surface is zero at these locations (x/c = 0.07 and 0.93) and varies smoothly in-between.
To ensure that the surface motion in the computational model is physically 210 realistic, a Bézier curve was used to define a C 1 continuous function from which the amplitude distribution was obtained. A Bézier curve is a parametric curve which is defined by the locations of 4 control points, 2 of which lie at each end of the curve, as seen in Figure 3. The frequency with which that surface moves is given by where η(x, t) is the local displacement in the direction normal to the undeformed 215 surface and V f + is the actuation frequency, non-dimensionalised by chord length, c, and freestream velocity, U ∞ .
(a) (b) Figure 4: Definition of the external boundary motion: a) Determination of amplitude as distance between points located along the normal to the cell face in the undeformed configurationn; b) Distribution of peak-to-peak amplitude along chord.
The temporal variation of the coordinates of the internal grid points is now considered. The purpose of the internal grid point motion is to accommodate the prescribed boundary motion as defined above, whilst preserving grid quality.
Recall that the grid velocity u pi appears in the equations of motion (1) and (2).
On the airfoil boundary u pi is equal to the prescribed velocity and at the external boundaries is equal to zero. A spatial distribution of u pi is therefore sought that provides acceptable grid distortion at each time instant. A smooth distribution is obtained if the grid velocities satisfy the Laplace equation, with a variable diffusion coefficient set equal to the square of the inverse of the distance from the moving surface, l, i.e.
The prescribed motion of the domain boundaries acts as a Dirichlet boundary condition for this equation.
The solution of (5) is computed iteratively prior to the external iterative 220 loop at each time-step. The Laplacian operator ensures that the displacement of the internal grid points is diffused smoothly inside the domain, as can be seen in Figure 5a. The grid quality in the near wall region must be maintained, so as little cell distortion as possible is desirable. As can be seen from Figure 5b, the maximum non-orthogonality is unaffected by the grid motion and the mean 225 non-orthogonality varies by less than 0.2% throughout the actuation cycle.

Results from static airfoil simulations
In this section the flow field around the static (i.e. unactuated) airfoil is characterised. Time-average results are first discussed, followed by analysis of the dominant flow structures obtained using the Dynamic Mode Decomposition 230 (DMD) method [33,34]. In all cases examined in this paper, the flow is homogenous in the spanwise direction, and the results are also spanwise-averaged to increase the sample size and improve convergence of statistics. When the flow is inhomogeneous in this direction (as in [35] due to spanwise periodic cambering, or in [36] due to spanwise flexibility) only time-averaging can be Upstream of x/c = 0.8 the levels of TKE are small, but downstream TKE increases rapidly, indicating turbulent transition. Horton [37] found that a reverseflow vortex and pressure recovery accompany the turbulent portion of a laminar separation bubble. Figure 6b shows that indeed a strong reverse flow vortex is confined to a region of high TKE. Moreover, a quiescent flow is expected to 250 develop beneath the laminar portion of the shear layer [38].  The flow development can be summarised as laminar separation with tran-sition but without reattachment. This is exactly the type of flow observed also experimentally [25,26]. There is also good quantitative matching between experiments and simulations. Velocity profiles at 3 locations in the airfoil wake are 275 shown in Figure 8.

Dominant structures in the separating shear layer
In order to understand the nature of the structures that govern the dynamics 285 of the shear layer in the uncontrolled flow field, DMD [33,34] was employed.
The method provides the most dynamically important structures, their spatial distribution and their frequency. For a thorough description of the method the reader is referred to Schmid [33], Jovanovic et al. [34] and Tu et al. [39].  The structures are still evident but they appear more corrugated due to the action of background turbulence.

Results from morphing airfoil simulations 310
In this section, the instantaneous, phase-averaged and time-averaged flow fields are analysed in order to provide insight into the flow control mechanism of surface morphing. As mentioned in Section 2.2, the actuation frequencies, V f , are non-dimensionalised as to give the reduced frequencies, V f + . Three frequencies were chosen for indepth investigation V f + = 0.4, 2 and 5, which correspond to V f = 14, 71 and 177 Hz respectively in the actual experiment [25,26]. The higher frequency was chosen to be close to the shear layer frequency identified in the previous section.
Due to structural constraints, this frequency could not be tested experimentally The flow was then allowed to develop for each V f + from the same starting conditions. The surface morphing assumes the form given by equation (4) in Section 2.2. The argument φ = 2πV f + t * in equation (6) is called the phase of 320 the actuation below. Note that because the deformation is fixed, the amplitude of the actuation velocity is proportional to the frequency V f + .  In both V f + = 2 ( Figure 10b) and V f + = 5 (Figure 10c), the first strong negative C f occurs at t * ≈ 0.8. Clearly the flow needs some time to respond to actuation. When V f + = 2, the regions are thicker and fewer than when Insight into the variation of the C f distributions is given by the instantaneous, spanwise-averaged Ω z contours shown in Figure 11. Figure  (a)   To investigate the structure of these vortices, Figure 14 displays instantaneous snapshots of the streamwise velocity, u/U ∞ , and an isosurface of Q-395 criterion at φ = 360 • for the 2nd and 4th cycles for V f + = 2 and 5 respectively.

Instantaneous Response of the Flow Field to Surface Motion
This time instant corresponds to the end of the actuation cycle in which shear layer roll-up first occurred, and it is equivalent to φ = 0 • for 3rd and 5th cycle in Figure 11. Q is defined as where Ω and S are the rotation and strain rate tensors respectively. When Q > 0 the 400 rotation rate dominates the strain rate, and serves to identify a vortex core. In As seen for the top plots of Figure 14, the flow fields still exhibit large regions of separated flow. As can be seen in Figure 15, this remains true even 4 cycles To investigate further the continued development of the flow field under the influence of control, the V f + = 2 simulation was continued until t * = 32.5 so it had also performed 25 cycles -equal to that of V f + = 5 case. The iso-contours of Ω z , u/U ∞ and the iso-surfaces of Q = 250 for φ = 360 • at the end of the 25th cycle are shown in Figure 16. Comparing the size of the separated regions in 420 Figure 16 with those in Figure 14 it is clear the actuation at V f + = 2 and 5 have significantly reduced the size of the separated region. When V f + = 5 the flow appears more spanwise coherent than when V f + = 2 and two LCS are present of the airfoil surface. The enhanced mixing that these produce has almost entirely eliminated the separated region at the trailing edge. Note also how close the 425 separating shear layer is in the airfoil (compare for example Ω z in Figure 16(b) and the plots in Figure 13).
To summarise, the control mechanism began with the early roll-up of the separated shear layer due to the amplification of the K-H instability. This first occurred during the 2nd and 4th cycle for the V f + = 2 and 5 respectively and  Figures 14-16, it has been found that this is a gradual process. As the number of actuation cycles increases, the shear layer gradually gets closer to the airfoil surface and the size of the separated region decreases.

Phase-Averaged Results
For a better understanding of the control mechanism, the spanwise-averaged Ω z and Q-criterion were phase-averaged every 45 • for the last 10 cycles. The temporal evolution of the flow fields is shown in Figures 18 and 19 for one complete cycle for V f + = 2 and 5 respectively.

450
For V f + = 2, the shear layer at the beginning of the cycle (Figure 18, φ = 0 • ) is located much closer to the surface than it was in Figure 11a, which is consistent with the large reduction in the size of the separation region observed in Figure 16. The Q-criterion reveals an LCS located at x/c ≈ 0.8. As the cycle continues to evolve, the LCS is ejected from the shear layer and propagates The fact that the spatial structure that forms at this frequency is similar 480 to that of the characteristic shear layer frequency (seen in Section 3.2) suggests that the periodic actuation is reinforcing a naturally occurring structure. This was a concept that was postulated by Benard and Moreau [19] who suggested that the actuation acts as a catalyser to amplify flow structures that already exist is the natural (i.e. unactuated) flow. Figure 18: Spanwise-averaged Ωz (Top) and Q-criterion (Bottom) phase-averaged over 10 actuation cycles when V f + = 2. Figure 19: Spanwise-averaged Ωz (Top) and Q-criterion (Bottom) phase-averaged over 10 actuation cycles when V f + = 5.
The phase-averaged plots reveal that once the flow has developed over the actuated airfoil, LCS are produced during every actuation cycle. This is consistent with the skin friction plots shown in Figure 10 and suggest that the shear layer is 'locked-on' to the surface motion. The lock-on of the shear layer to the forcing frequency was also observed by Kotapati et al. [13] in their investigation 490 of synthetic-jet based periodic control.
The following mechanism is then proposed: the triggering of K-H instability is the first in a chain of events that eventually leads to separation control.
Amplification of the K-H instability causes the early onset of shear layer rollup and the creation of LCS. By entraining fluid from the outer flow and re-to the airfoil surface with each actuation cycle. It takes at least 15 cycles for this process to complete. Once this happens, the structures are directly re-energising the near-wall flow, that recovers enough energy to overcome the adverse pressure gradient. Consequently the separation zone is significantly 500 reduced. This mechanism will be further explored in the following section.

Time-Averaged Results
The process of time-averaging a velocity field involves averaging the veloc-   For the frequencies of V f + = 2 and 5, the time-averaging was performed us-525 ing 80 snapshots from the final 10 cycles (i.e. 8 snapshots per cycle). Figure 21 demonstrates that this number of samples is sufficient to provide a statistically converged time average for both actuation frequencies. The time-averaged streamwise velocity, TKE, and Reynolds stress, u v , are shown in Figure 22.
(c)V f + = 5 In Section 3.1 it was found that for the static airfoil, transition to turbulence 530 occurred as a result of the natural breakdown of the separated shear layer. This means that the regions of high TKE and Reynolds stress occur in the shear layer, away from the airfoil surface. This can be seen in Figure 22 can therefore be confirmed that dynamic surface motion at V f + = 2 and 5 is 545 able to gradually move the region of high momentum transfer close to the airfoil surface. This then results in significant reduction in separation.

Effect of Actuation on Aerodynamic Performance
The effect that actuation on the aerodynamic performance of the airfoil can be seen from the C P distribution in Figure 23  As can be seen in Table 1, this results in improvements in the aerodynamic force coefficients. When V f + = 2, C L increases by 183% and C D decreases by 560 32%. V f + = 5 has a greater decrease in C D with a 37% reduction but the lift recovery is slightly less at 141%. This equates to a threefold L/D increase from 2.8 to 9.4 and 8.8 for V f + = 2 and 5 respectively. In both frequencies the flow is locked to the wall oscillation. Similar lift enhancement has been also observed when the motion of flexible membrane wings is synchronised with 565 vortex shedding [44,45].
The C P plot indicates that frequency affects not only the pressure distribution on the suction side, but also the distribution on the pressure side. The combined effect results in the small drop in C L and the ratio C L /C D when V f + increases from 2 to 5. It is not immediately evident how the top wall oscillation 570 affects the pressure at the bottom wall. The effect is probably inviscid, because the flow is still attached in this area. It may be related to the periodic acceleration and deceleration of the mass adjacent to the moving wall, that has a global effect due to pressure coupling (continuity equation). The morphing wall covers 86% of the suction side, so a large mass is displaced. Interestingly, in the paper 575 of Kang et al [24] the same behaviour is also observed, when a small part of the wall in the suction side (close to the leading edge) is oscillating. In the lock-in region, the lift does indeed increase, but the increase is not monotonous and the authors do not report the ratio C L /C D . The pressure at the bottom wall is also affected by the motion.   the viscous drag is substantially increased. Often the decrease in the C DP is greater than the associated increase in C Dv and therefore a net reduction in C D does occur. The particular control method examined in this paper does not rely on tripping the boundary layer. From Table 1 it can be seen that the contribution of pressure drag decreases and that of viscous drag slightly increases.

590
The increase in the viscous component is a consequence of the increase in skin friction from the stronger vortex shedding. In absolute terms however the pressure drag, C Dp , is reduced by 20% and 31% for V f + = 2 and 5 respectively compared to the static case. There is an associated relatively small increase in viscous drag.

Conclusions
In this paper the effect of dynamic surface motion on the flow development around the airfoil was explored computationally. The flow separated from the airfoil surface, transitioned but did not reattach. Actuation frequencies were chosen with reference to the fundamental shear layer frequency for the unac-600 tuated airfoil that originates from Kelvin-Helmholtz instability. It was found that after a few actuation cycles the Kelvin-Helmholtz instability was amplified and the shear layer started to roll-up earlier. LCS were formed in the shear layer which broke down and enhanced momentum transfer from the outer flow.
As these were repeatedly formed and propagated downstream, the enhanced and Reynolds stresses revealed that this resulted in the region of high momentum exchange occurring in the near-wall region, thus re-energising directly the boundary layer. The overall result was an almost complete elimination of the separated regions and a significant improvement in aerodynamic performance.

615
The first author acknowledges financial support from the International Joint-PhD Scholarship of the International Office of Imperial College. Significant financial assistance was also provided by Temasek Laboratories, National University of Singapore, Singapore. The simulations were carried out on the facilities of Archer, the UK national high-performance computing service, under the 620 grant EP/L000261/1 of the UK Turbulence Consortium (UKTC) and on the High Performance Computing Service at Imperial College London (CX2).