Nomenclature

c

chord length

C_t

thrust coefficient

C_P

power coefficient

dy

wall distance

f

plunge frequency

h

plunge position

h ₀

plunge amplitude

k

reduced frequency, $= \omega {c}/{U}_{\infty}$

Re

Reynolds number, $\rho {U}_{\infty} {c}/\mu$

Re ${}_{\theta}$

momentum thickness Reynolds number

${U}_{\infty}$

freestream velocity

v_p

non-dimensional plunge velocity, $={k}({h}_0/{c})$

ITA

iterative time advancement scheme

MAV

micro-air-vehicle

OVG

oscillating vortex generators

PISO

pressure-implicit with splitting of operator

SST

shear-stress transport

URANS

unsteady Reynolds-averaged Navier-Stokes

VG

vortex generators

$\gamma$

intermittency

$\eta$

propulsive efficiency, $= C_{{t}}/C_{{P}}$

$\rho$

fluid density

$\omega$

angular frequency, $=2{\pi f}$

2.0 Computational setup

Dohring [Reference Dohring15] used a 760mm long flat plate aligned with the water flow and a NACA 0012 aerofoil of 20mm chord length c whose leading edge was located 100mm downstream of the plate leading edge, Fig. 2(a). The computations are performed as the aerofoil undergoes a periodic plunging motion given by Equation (1) The geometric angle-of-attack of the aerofoil is set to 0°. The distance between the flat plate and the plunging NACA0012 aerofoil is called the wall distance (dy). It ranges between 0.4c to 1.5c, corresponding to 8 to 30mm. The plunging amplitude h of the aerofoil is 0.088c, corresponding to 1.76mm.

Figure 2(b) shows the schematic diagram of counter-rotating biplane aerofoil near the flat plate. This configuration is built to illustrate the MAV shown in Fig. 1. The two aerofoils are again flapping in counterphase, but they are mounted upstream of the plate trailing edge. The aerofoil distance from the plate was varied up to a distance where the foils were flapping outside of the plate boundary layer. The flat plate has semi-circular rounded leading and trailing edges, and its thickness and length are set to 0.5c and 5c, respectively.

The sinusoidal plunging position (h) of the aerofoil and the angular frequency ( $\omega$ ) are defined by Equations (1) and (2), respectively.

(1) \begin{align}h\left( t \right) &= {h_0}sin\left( {\omega t} \right)\end{align}

(2) \begin{align}\omega &= 2\pi f\end{align}

where h ₀ is the plunge amplitude, and f is the plunge frequency in Hz. Other critical parameters for the plunge motion are the reduced frequency (k) (Equation 3) and non-dimensional plunge velocity ( $v_{p}$ ) (Equation 4).

(3) \begin{align}k = \omega c/{\rm {U_\infty}}\end{align}

(4) \begin{align}{v_p} = k\left( {{h_0}/c} \right)\end{align}

where c is the chord length, and U _∞ the freestream velocity. The propulsive efficiency ( $\eta $ ) is then defined by the ratio between the mean thrust (C _t ) and mean power (C _P ) coefficients.

(5) \begin{align}\eta = {C_t}/{C_P}\end{align}

The mean power coefficient is computed by integrating the lift coefficient times the plunge velocity on the plunging aerofoil, as shown in Equation (6). Time-averaged data is obtained by taking the mean of 10 complete oscillations after periodic conditions are reached.

(6) \begin{align}{{{C}}_{{P}}} = \frac{1}{{\Delta t}}\mathop \smallint \limits_{\Delta t}^{\rm{\;}} \left( {{c_l}\dot h/{U_\infty }} \right){\rm{d}}t\end{align}

Figure 2. Schematic diagram of (a) single plunging aerofoil near the flat plate and (b) counter-plunging biplane aerofoil near the flat plate.

3.0 Numerical method

The two-dimensional low-Reynolds-number incompressible flow over the plunging NACA0012 near the flat plate is simulated numerically with a commercial finite volume solver of ANSYS Fluent v17.2 [17] by solving the unsteady Reynolds-averaged Navier-Stokes (URANS) equations. The numerical solution is obtained by evaluating the gradients using the Least Squares Cells-Based scheme. The second-order upwind scheme was adopted as the method to discretise pressure, momentum, turbulent kinetic energy, specific dissipation rate, intermittency and momentum thickness Re equations. The transient formulation uses the second-order implicit iterative time advancement scheme (ITA). The advantage of the implicit scheme is that it is unconditionally stable with respect to time-step size [17]. Pressure-velocity coupling is done by a segregated type algorithm of PISO (pressure-implicit with splitting of operators) and the details of the algorithm can be found in Ferziger et al. [Reference Ferziger, Peric and Street18]. For the pure plunging simulations of a single NACA0012 aerofoil, the motion of the aerofoil is introduced by rigid-body motion of the whole domain. The high-quality structured grid near the NACA0012 aerofoil (shown in Fig. 3 as red region named inner rigid zone) moves together with the aerofoil as a rigid body, and the blue-coloured region is a deformable grid. Other regions coloured as grey are set as stationary and non-deformable (see Fig. 2).

Figure 3. Computational domain.

Figure 4. Time history of the thrust coefficient at (a) Re = 1,960, k = 12.7, dy = 0.4c and (b) Re = 5,940, k = 4.2, dy = 0.4c.

The numerical simulations were conducted with the transitional turbulence model of the shear-stress transport (SST) $\gamma $ - $R{e_\theta }$ model (also known as Transition SST model), which couples the $k$ - $\omega $ SST model of Menter [Reference Menter19] with two transport equations for intermittency ( $\gamma $ ) and momentum thickness Reynolds number ( $R{e_\theta }$ ) (Langtry and Menter [Reference Langtry and Menter20]). The model has been used and validated at low Reynolds number external flow over an aerofoil in the literature by Counsil and Boulama [Reference Counsil and Boulama21], Chen et al. [Reference Chen, Yu and Mei22] and Ruiz and D’Ambrosio [Reference Ruiz and D’Ambrosio23]. Time-averaged (mean) data is obtained by taking the mean of a sufficiently large number of complete oscillations after periodic conditions are reached.

The velocity inlet is placed sufficiently far upstream to ensure uniform freestream flow at the beginning of the flat plate. The inlet is 1 $c$ upstream of the leading edge of the flat plate. The bold blue lines at the top and bottom boundaries show the symmetry boundary conditions. Also, the bold red line (total 38 $c$ length) at the bottom boundary indicates the flat plate with wall boundary condition.

3.1 Grid refinement study

In our previous papers [Reference Anilir, Kurtulus and Platzer16,Reference Anilir, Kurtulus and Platzer24], the grid refinement study was carried out by adjusting the number of points around a plunging NACA0012 aerofoil and the structured cells at the near wake. Therefore, the inner rigid zone (red region) shown in the Fig. 3 is the same as in the previous studies and the number of points around the aerofoil surface is equal to 760. For this study, the grid refinement was only applied to the cell on the flat plate surface at Re = 1,960, k = 12.7, dy = 0.4c, and Re = 5940, k = 4.2, dy = 0.4c. Three different grids, named as coarse, medium and fine grids, were employed with a refinement ratio of around 1.5. The time history of the plunging aerofoil’s thrust coefficient is presented in Fig. 4. It is found that the medium grid is sufficiently refined and the results are very close to the fine grid. The step size of 0.0002s (500 time steps within a plunge cycle) based on the previous time refinement studies is selected.

Figure 5. Comparison of the mean streamwise velocity profiles at three wake distances for k = 57, v_p = 5.03, Re = 440, dy = 0.4c.

4.0 Validation of the numerical simulations

The numerical simulations are validated against the published experimental results of Dohring [Reference Dohring15] at several reduced frequencies, Reynolds numbers and wall distances. The experiments were carried out in a water tunnel by conducting Laser Doppler measurements and flow visualisation for the case of a plunging aerofoil in close proximity to a flat plate. Figure 5 shows the comparison between the computed and measured velocity distributions at three different measuring stations downstream of the foil trailing edge for a freestream velocity of 0.022m/s, corresponding to a Reynolds number of 440 based on the foil chord. The laminar flat plate boundary layer thickness at the foil leading edge (located 100mm downstream from the leading edge of the flat plate) then is 10mm. The foil distance from the flat plate, the wall distance dy, is 8.5mm, corresponding to a non-dimensional wall distance dy = 0.425c. The aerofoil therefore is just barely immersed in the boundary layer. The amplitude of oscillation is 1.76mm and the frequency of 10Hz corresponds to a reduced frequency of 57. The time-averaged velocity profiles downstream of the aerofoil trailing edge are shown at distances 135, 150 and 200mm from the plate leading edge, corresponding to local boundary layer thicknesses of 12.4, 13 and 16.9mm, respectively. The formation of a jet profile is clearly seen, decreasing from a value of more than eight times the freestream velocity at a distance of 1.75c to four times at a distance of 7.5c. These profiles are predicted quite well. At large downstream distances the jet has dissipated into the freestream.

In Fig. 6, the comparisons are shown for two much larger wall distances dy = 1.0c and 1.5c, corresponding to distances of 20 and 30mm, respectively. Consequently, the aerofoil is outside of the boundary layer. The free-stream speed is 3m/s corresponding to a Reynolds number of 5,940, the frequency is 10Hz corresponding to a reduced frequency k = 4.2. As expected, the jet profiles are much smaller than in Fig. 5. The agreement between experiment and computation is quite good.

Figure 6. Comparison of the mean streamwise velocity profiles at the wall distances dy = 1.0 and dy = 1.5 for k = 4.2, v_p = 0.37, Re = 5,940.

Similarly, the velocity profiles are shown in Fig. 7 for the wall distances of dy = 0.4c and 1.25c, corresponding to distances of 8 and 25mm, respectively. The free stream speed is 0.1m/s corresponding to a Reynolds number of 1960 and the frequency is 10Hz corresponding to a reduced frequency of 12.7. It is seen that the maximum velocity is larger at the larger wall distance compared to the maximum velocity at the much smaller wall distance. This indicates that there exists a maximum velocity for different parameter combinations of wall distance and plunge velocity. A comparison of the optimal wall distances obtained from the measurements and computations is given in Table 1. It is seen that the optimal distance increases with increasing plunge velocity. The comparison is shown in Fig 8 for four plunge velocities. It is quite good with the exception of a few points in Fig. 8(b), (c) and (d) that require further investigation.

Furthermore, it is also interesting to note the upstream effect in front of the foil leading edge, shown in Fig. 9 At the two distances 0.4c and 0.75c in front of the leading edge, corresponding to 8 and 15mm, the flow speed increases in both the experiment and computation. The differences are attributable to the measuring inaccuracy.

Finally, Fig. 10 shows the comparison between Dohring’s flow visualisation of the vortical wakes shed from the trailing edge of the plunging foil for a plunge velocity of 0.37 at the two wall distances dy = 0.6c and 1.1c. The agreement is quite good. At the closer distance the vortical wake is deflected upward from the wall. At the larger distance the vortical wake remains symmetrical in agreement with the expected behaviour in unbounded flow at this plunge velocity.

Table 1. Optimal wall distances

Figure 7. Comparison of the mean streamwise velocity profiles at the wall distances dy = 0.4c and 1.25c for k = 12.7, $v_{p}$ = 1.12, Re = 1,960.

Figure 8. Variation of the nondimensional maximum velocity in the velocity profiles obtained at x = 1.75c at (a) k = 4.2, (b) k = 6.3, (c) k = 8.4, and (d) k = 12.7.

Figure 9. Comparison of measured and computed velocity profiles upstream of the foil leading edge for k = 57, dy = 0.6c.

Figure 10. Comparion of vortical patterns between numerical simulations and Dohring’s water tunnel experiments $v_{p}$ = 0.37, (a) dy = 0.6c, and (b) dy = 1.1c.

Figure 11. Variation of thrust coefficient, lift coefficient and propulsive efficiency for (a) k = 12.7, $v_{p}$ = 1.12, Re = 1,960, (b) k = 8.4, $v_{p}$ = 0.74, Re = 2,990, (c) k = 6.3, $v_{p}$ = 0.55, Re = 3,960, and (d) k = 4.2, v_p = 0.37, Re = 5,940.

5.0 Computational prediction of the aerodynamic characteristics of plunging naca 0012 profiles as a function of wall proximity

For the application of small oscillating aerofoils to boundary layer propulsion, it is necessary to know their thrust, lift and efficiency. In this section the above described URANS simulations are used to predict the aerodynamic characteristics of NACA 0012 profiles oscillating in the pure plunge mode at various distances from the wall. Figure 11 shows the variation of the thrust and lift coefficients and of the propulsive efficiency for wall distances between 0.5c to 5.0c at various plunge velocities and Reynolds numbers. The dashed line indicates the values for the foil in unbounded flow. It is seen that the thrust force produced by the plunging aerofoil near the flat plate is always greater than that of the single plunging aerofoil in unbounded flow, but continuously decreases with increasing distance to approach the unbounded flow value. A similar behaviour is observed for the lift and the propulsive efficiency. As expected, the thrust increases with increasing plunge velocity, whereas the propulsive efficiency exhibits a consistent peak at a plunge velocity of about 0.55, as shown in Fig. 12.

Figure 12. Variation of (a) the thrust coefficient and (b) propulsive efficiency with respect to nondimensional plunge velocity for different wall distances.

Figure 13. Instantaneous vorticity field for (a) dy = 0.4c, (b) dy = 0.5c, (c) dy = 0.6c, (d) dy = 0.75c, (e) dy = 0.9c, (f) dy = 1.0c, (g) dy = 1.1c, (h) dy = 1.5c, and (i) single plunging aerofoil at k = 8.38, $v_{p}$ = 0.74, Re = 2,990.

It is instructive to observe the change in vortex shedding as a function of wall distance and plunge velocity. Dohring [Reference Dohring15] found it necessary to distinguish between wall distances dy<0.5c, 0.5c<dy<1.0c, dy>1.0c and plunge velocities $v_{p}$ >1 and $v_{p}$ <1. It is well known that in unbounded flow the vortex shedding occurs in the form of a symmetric reverse Karman vortex street for $v_{p}$ <1 that becomes asymmetric for $v_{p}$ >1, deflected upwards or downwards depending on the initial starting motion of the foil. Dohring’s flow visualisation revealed that for dy<0.5c and small $v_{p}$ a row of clockwise vortices is shed along the plate and counterclockwise vortices are weak or not visible at all. For $v_{p}$ >1 a row of counterclockwise vortices is shed along the plate and the clockwise vortices are weak or invisible. Therefore, in the former case the induced velocities are negative and the maximum flow velocity cannot occur near the plate, whereas in the latter case the induced velocities are positive and the maximum velocity will always be close to the plate. For distances between 0.5c and 1.0c, the ground effect changes the symmetric vortex shedding into an upward deflected shedding of vortex pairs. This change is visible in Fig. 13(a), (b), (c), (d) for wall distances < 0.9c. The deflection does not depend on the initial plunge direction. At wall distances >1.0c, its effect has become minimal. This pairing phenomenon has also been observed for pure pitch oscillations in ground proximity by Quin et al. [Reference Quinn, Moored, Dewey and Smits25]. For low k and small wall distances, the vorticity is elevated away from the ground and pushed downstream by the counterclockwise vortex as shown in Fig. 14. This vorticity is denoted as ground vorticity by Molina et al. [Reference Molina, Zhang and Angland26].

6.0 Application to flat plate boundary layer flow

The MAV, shown in Fig. 1, benefited greatly from the ground effect phenomenon by positioning two aerofoils flapping in counterphase in a biplane arrangement downstream of the fixed wing’s trailing edge. In order to demonstrate the direct boundary layer propulsion method by oscillating foils proposed in the present paper, the flow over the flat plate configuration shown in Fig. 2(b) was computed. The two foils are again flapping in counterphase, but they are mounted upstream of the plate trailing edge.

Figure 15 shows the variation of the thrust produced by the upper aerofoil, lower aerofoil, the drag of the flat plate and the thrust of the complete configuration at the wall distances of 0.4c and 1.5c by adding the thrust of the two oscillating foils and subtracting the drag of the plate. The coefficients are nondimensionalised with the chord length of the aerofoils. For k = 12.7, the upper and lower aerofoils generate an almost identical thrust force. As expected, their values decrease from around 1.1 to 0.75 as the wall distance increases from 0.4c to 1.5c. These values are very close to the results obtained in the previous section from the plunging aerofoil near the flat plate at a position far upstream from the plate trailing edge. However, a significant difference can be observed in Fig. 16 between oscillating the foils in the boundary layer at dy = 0.4c versus outside the boundary layer at dy = 1.5c. As expected, no interaction takes place between the Karman vortex street shed from the flat plate in the latter case. In the former case, the Karman vortex street is prevented from forming, resulting in an approximately 40% drag reduction of the flat plate.

Similarly, as shown in Fig. 17, for the lower plunge frequency and higher Reynolds number case of k = 4.2 and Re = 5,940, there is no interaction between the Karman vortex street and the vortices shed from the aerofoils at a wall distance of dy = 1.5c, but a very strong interaction at dy = 0.4c. In the former case the thrust of the aerofoils is too small to overcome the plate drag while in the latter case it is still sufficient to produce a net thrust.

Figure 14. Instantaneous vorticity field for dy = 0.6c at (a) k = 6.3, and (b) k = 4.2.

Figure 15. Influence of wall distance on the thrust produced by upper aerofoil, lower aerofoil, and flat plate at (a) k = 12.7, and (b) k = 4.2.

Figure 16. Instantaneous vorticity field for biplane arrangement for k = 12.7, $v_{p}$ = 1.12, Re = 1,960 at (a) dy = 0.4c, and (b) dy = 1.5c.

Figure 17. Instantaneous vorticity field for biplane arrangement for k = 4.2, $v_{p}$ = 0.37, Re = 5,940 at (a) dy = 0.4c, and (b) dy = 1.5c.

7.0 Summary and outlook

The flight speed of vehicles using flapping wing propulsion is limited because the required high Strouhal number can only be achieved by a low speed in the denominator of the Strouhal number. However, aerofoils or hydrofoils oscillating in a vehicle’s boundary layer encounter low flow speeds while the vehicle’s flight speed may be significantly larger. The investigation of oscillating foils in or close to vehicle boundary layers therefore may open possibilities for boundary layer energisation that deserve exploration.

The analysis presented in this paper shows that small foils oscillating in the pure plunge mode in or close to the boundary layer of a flat plate generate significantly larger thrust than foils oscillating in an unbounded flow of equal Reynolds number, thus indicating the possibility of favorably influencing the boundary layer behaviour on large-scale aerofoils or hydrofoils. The results also show beneficial interaction effects between the Karman vortex street shed by the flat plate and the vortices shed by the flapping foils. Such foils therefore can be regarded as oscillating vortex generators (OVG) in generalisation of the well-known steady-state (VG) vortex generators.