Universal scaling law in drag-to-thrust wake transition of flapping foils

TL;DR

This paper introduces a universal scaling law for the wake transition from drag to thrust in flapping foils, based on the cycle-averaged swept trajectory, providing a unified framework for understanding propulsion mechanisms.

Contribution

The study replaces amplitude with swept trajectory in the scaling law, revealing a universal $St_{\mathcal{T}} \sim 1$ criterion for wake transition across different flapping modes.

Findings

Wake transition occurs at $St_{\mathcal{T}} \sim 1$ for all tested modes.

Propulsive jets form when the induced velocity exceeds free stream velocity.

The new metric links wake development directly to kinematic parameters.

Abstract

Reversed von K\'arm\'an streets are responsible for a velocity surplus in the wake of flapping foils, indicating the onset of thrust generation. However, the wake pattern cannot be predicted based solely on the flapping peak-to-peak amplitude $A$ and frequency $f$ because the transition also depends sensitively on other details of the kinematics. In this work we replace $A$ with the cycle-averaged swept trajectory $\mathcal{T}$ of the foil chord-line. Two dimensional simulations are performed for pure heave, pure pitch and a variety of heave-to-pitch coupling. In a phase space of dimensionless $\mathcal{T}-f$ we show that the drag-to-thrust wake transition of all tested modes occurs for a modified Strouhal $St_{\mathcal{T}}\sim 1$. Physically the product $\mathcal{T}\cdot f$ expresses the induced velocity of the foil and indicates that propulsive jets occur when this velocity exceeds $U_{\infty}$. The new metric offers a unique insight into the thrust producing strategies of biological swimmers and flyers alike as it directly connects the wake development to the chosen kinematics enabling a self similar characterisation of flapping foil propulsion.

Figures

Tables

Table 1. Table 1: Computational statistics of grid convergence for a harmonically flapping foil.

Grid density	Thrust coefficient	$\left|\Delta C_t\right|$	Relative $\left|\Delta C_t\right|$ %
064	0.494	9.3E -2	15.8
128	0.555	3.2E -2	5.45
192	0.575	1.2E -2	2.04
256	0.580	0.7E -2	1.19
512	0.587	—	—

Equations

$$\begin{array}[]{l}\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}y_{c}(t)=y_{h}(t)+y_{\theta}(t)~{}~{}~{}where:\\ \\ \displaystyle y_{h}(t)=h_{\mathit{0}}\sin(2f\pi t)~{}~{}~{}~{},~{}~{}~{}~{}y_{\theta}(t)=(1-\mathcal{P})\mathcal{C}\sin({\theta}(t))~{}~{}~{}and:\\ \\ \displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\theta(t)=\theta_{\mathit{0}}\sin(2f\pi t+\psi)\\ \end{array}$$

$$\alpha =\arctan \frac{\dot{h_0}}{U_{\infty }}-{\theta }_0$$

$$Sr=\frac{Df}{U_{\infty }},A_D=\frac{2y_t(t_{max})}{D}$$

$$C_t=\frac{F_x}{\frac{1}{2}\rho U_{\infty }^{2}c}$$

$$T_r=U_{\infty }{\int }_0^{1/f}\sqrt{1+{[\frac{d{y}_t(\mathcal{C},t)}{dt}]}^2}dt$$

$$\mathcal{T}=\overline{Tr}=U_{\infty }{\int }_0^1{\int }_0^{1/f}\sqrt{1+{[\frac{d{y}_t(s,t)}{dt}]}^2}dtds$$

Full text

1 Introduction

Almost all aquatic and flying animals generate thrust via the oscillatory motion of foil-like body parts e.g. tail, fin etc. Moreover flapping foil systems are often associated with high efficiency and strong side forces, ideal for manoeuvring (Read et al., 2003). Thus, many studies have focused on the analysis and implementation of these biological configurations into man made designs (Fish & Lauder, 2006; Triantafyllou et al., 2004; Wang, 2005) although the underlying physics is still not clearly understood.

Here we aim to determine the drag-to-thrust wake transition of these flapping mechanisms via the analysis of the vortex pattern development. Assuming foil undulations above the lock in frequency (Thiria et al., 2006; Vial et al., 2004) we observe at least three basic wake patterns (Von Karman, 1935): the classic $B\acute{e}nard~{}\emph{von}~{}K\acute{a}rm\acute{a}n$ (Bvk) street where $U_{wake}<U_{\infty}$ (figure 1a), the neutral line where $U_{wake}\sim U_{\infty}$ (figure 1b) and the reversed BvK wake, where $U_{wake}>U_{\infty}$ (figure 1c). The latter is synonymous to the drag-to-thrust wake transition although a lag exists between this phenomenon and the foil’s overall transition towards thrust. This is due to the fact that a weak velocity surplus cannot overcome profile drag or velocity fluctuations and pressure differences within the control volume(Streitlien & Triantafyllou, 1998; Ramamurti & Sandberg, 2001; Bohl & Koochesfahani, 2009).

As the driving factors of BvK reversal we typically consider the oscillating amplitude and the oscillating frequency $f$ of the kinematics (Koochesfahani, 1989) . The former is expressed by the trailing-edge (TE) peak-to-peak amplitude $A$. In dimensionless terms it is often normalised by the thickness D or the chord length $\mathcal{C}$ of the foil. In a similar fashion the frequency is often expressed as a reduced frequency $k=U_{\infty}/(f\mathcal{C})$ (Birnbaum, 1924), a thickness based Strouhal number $Sr=(fD)/U_{\infty}$ (Godoy-Diana et al., 2009) or a chord length based Strouhal number $St_{\mathcal{C}}=1/k$ (Cleaver et al., 2012). Triantafyllou et al. (1991) suggested a modified amplitude based Strouhal number $St_{A}=(fA)/U_{\infty}$. By including both the frequency and the amplitude of oscillation, $St_{A}$ can potentially characterise the BvK reversal by a single factor as opposed to $k$, $Sr$ and $St_{\mathcal{C}}$. Studies of Anderson et al. (1998) and Read et al. (2003) showed that optimal efficiency occurs for a short range of $St_{A}\sim[0.2,0.4]$. This was also supported by Taylor et al. (2003) and Triantafyllou et al. (1993) who observed that the majority of natural fliers and swimmers prefer to cruise within this range. According to Andersen et al. (2017), BvK reversal occurs at different $St_{A}$ values for pure heave and pure pitch. Therefore, $St_{A}$ cannot be regarded as an expression of self similarity.

The fundamental problem is that characterizing the motion only by the tail amplitude fails to capture the contribution of the other points of the foil. This becomes important when the heaving component is significant, resulting in the generation of strong leading edge vortices (LEV) which travel downstream and blend with the trailing edge vortices (TEV). Instead, we need to take into account the length of the entire path travelled in a period rather than the maximum distance from equilibrium expressed by $A_{D}$.

In this study, we formulate a novel length scale, which characterizes the Bvk reversal of harmonically flapping foils. Two-dimensional simulations are conducted at a Reynolds number of $Re=1173$ for a rigid NACA0016 profile and three basic harmonic kinematics: pure heave, pure pitch and heave-pitch coupling. Additional higher Reynolds number simulations ($Re=11730$) are used to quantify the Reynolds number effects. The influence of different pivot points is examined for pure pitch and coupled motions. In addition, we analyse the impact of different maximum effective angles of attack on coupled kinematics. Finally, we develop a new metric based on the chord-wise averaged path travelled by the foil, in order to determine BvK wake reversal for a vast range of harmonic motions.

1.1 Geometry and Kinematics

We consider a rigid NACA0016 with a thickness $D=0.16\mathcal{C}$. The foil performs simple harmonic oscillations around a stationary equilibrium position, against a uniform free stream velocity $U_{x}=U_{\infty}$. The lateral direction of every point along this chordline is denoted with $y(t,s)$. Here $t$ is the time and $s$ is the $\mathcal{C}$-normalised coordinate along the chord ranging from 0 at LE to 1 at TE. Pure pitch is modelled as a sinusoidal rotation about a specified pivot point along the chordline ($s=\mathcal{P}$, where $\mathcal{P}$ is the non-dimensional distance between the LE and the pivot point along the chordline) and pure heave as a sinusoidal lateral translation of the entire chordline. Coupled motion occurs by the superposition of these pure motions:

[TABLE]

Here $\theta(t)$ is the instantaneous value of pure pitch whilst $h_{0}$ and $\theta_{0}$ are the amplitudes of pure heave and pure pitch respectively. The phase difference between pitch and heave is expressed as $\psi$. Here $\psi$ = $90^{\circ}$ because this value is considered optimal in terms of propulsive efficiency (Platzer & Jones, 2008).

Another important kinematic parameter in coupled motions is the effective angle of attack $\alpha_{\emph{eff}}(t)$ which is the summation of the instantaneous pitch angle $\theta(t)$ and the heave induced angle of attack. Thus for $\psi$ = $90^{\circ}$ the amplitude of $\alpha_{\emph{eff}}(t)$ is:

[TABLE]

where $\dot{h_{0}}$ is the amplitude of $\dot{y_{h}(t)}$ = $dy_{h}/dt$. In this study we differentiate coupled motions by varying $\mathcal{P}$ and $\alpha$ while keeping $\psi$ = $90^{\circ}$.

1.2 Dimensionless parameters

Three non-dimensional parameters are used to describe the interaction between the oscillating foil and the free stream: the Reynolds number based on chord length, $Re=U_{\infty}\mathcal{C}/\nu$ (where $\nu$ is the kinematic viscosity) the thickness based Strouhal number, $Sr$ and the non-dimensional TE amplitude ($A_{D}$):

[TABLE]

Thus $St_{A}=Sr\cdot A_{D}$ and can be understood as the ratio between the speed of the foil tip and $U_{\infty}$ (Godoy-Diana et al., 2009).

The thrust coefficient $C_{t}$ is expressed by simply normalising the total force acting on the $x-\emph{axis}$ by the dynamic pressure of the freestream. Time averaged quantities are presented with an overbar to distinguish them from their instantaneous counterparts.

[TABLE]

1.3 Computational Method

The CFD solver chosen for this study can simulate complex geometries and moving boundaries for a wide range of Reynolds numbers in 2D and 3D domains, by utilizing the boundary data immersion method BDIM, (Weymouth & Yue, 2011). BDIM solves the viscous time-dependent Navier-Stokes equations and simulates the entire domain by combining the moving body and the ambient fluid through a kernel function. This technique has quadratic convergence and has been validated for flapping foil simulations over a wide range of kinematics (Maertens & Weymouth, 2015; Polet et al., 2015).

The mesh configuration is a rectangular Cartesian grid with a dense uniform grid near the body and in the near wake, and exponential grid stretching used in the far-field and the numerical domain uses a uniform inflow, zero-gradient outflow and free-slip conditions on the upper and lower boundaries.Furthermore no slip boundary conditions are used on the oscillating foil. Mesh density is expressed in terms of grid points per chord. A uniform grid of $\delta x=\delta y=\mathcal{C}/192$ is used for the results in this work based on the results of the convergence study shown in table 1.

2 Results and Discussion

2.1 Wake comparison of different kinematics

According to Godoy-Diana et al. (2008), the Reynolds number range of naturally occurring flapping foils is $100<Re<10000$. In this study, most simulations are conducted for $Re=1173$ to be within this range and to enable comparison with the experiments of Godoy-Diana et al. (2009). Additional simulations are conducted for selected cases at $Re=11730$ to examine the higher Reynolds number effects. The pivot points tested for pure pitch and coupling are $\mathcal{P}=0$ and $0.25$. The coupled motion is also tested for three values of $\alpha=5^{\circ},10^{\circ}$ and $20^{\circ}$.

We analyse the wake patterns and resultant hydrodynamic loads for the above mentioned kinematics. Various stages of the wake development can be seen in figure 3 for pure heave, pure pitch for $\mathcal{P}=0.25$ and coupled motion for $\mathcal{P}=0.25$ and $\alpha=10^{\circ}$. The transition from the BvK (third row) to the neutral wake where vortices are shed in-line (fourth row) and later the reversed BvK (fifth and sixth rows) is in agreement with literature (Koochesfahani, 1989; Godoy-Diana et al., 2009; Andersen et al., 2017) . At lower $Sr-A_{D}$ combinations more complicated wake patterns e.g. 2P wakes (Williamson & Roshko, 1988) at the first row of figure 3 in accordance with those observed by Andersen et al. (2017) for wedge type foils. At such low $Sr-A_{D}$ combinations coupled motions are dominated by one of the two modes e.g. for $A_{D}=0.4,Sr=0.1$, the heaving contribution to the foil displacement is less than $5\%$ for all the coupled cases studied. Therefore, a coupled motion around this region is acting more like a pure pitching case.

Among the three kinematic test cases significant discrepancies can be seen in close proximity to the foil, most notably at the LE. The deep stall (high $\delta\alpha/\delta t$ across the chord) experienced by the pure heaving foil generates LEVs of sizes comparable to the TEVs which travel across the chord and blend with the wake. A closer look at figure 3c reveals that the coupled motion generates the smallest amount of dynamic separation among the three cases. Finally as seen in the last row of figure 3 even when the BvK street is fully reversed, some the cases exist within the drag producing regime. This lag is expected since a weak propulsive wake is not enough to overcome the profile drag or to compensate for the velocity fluctuations and pressure differences that exist within the control volume (Streitlien & Triantafyllou, 1998; Ramamurti & Sandberg, 2001; Bohl & Koochesfahani, 2009).

Figure 4 shows the best fit curve that isolates the neutral line (where $U_{wake}\sim U_{\infty}$), for the different harmonics. These best fit curves are reproduced in figure 5a in order to compare this neutral line across different kinematics. Although $Sr-A_{D}$ phase diagram is suitable to examine a specific kinematics, it is clear that this does not universally describe wake transitions. This is the result of the unique interactions between LEVs and TEVs that are specific to the motion type. Consequently this demonstrates the need for a self similar classification of the oscillating amplitude to accurately determine wake development.

2.2 An Alternative Length Scale

Fundamentally, the foil generates thrust force by displacing and accelerating fluid out of its path as it moves through its prescribed trajectory. The quantity of fluid displaced is dependent on the product of the chord length $\mathcal{C}$ times the path length travelled by the foil over one period of oscillation. Thus,a proper indicator of the wake’s drag-to-thrust transition should reflect the length of the curve traversed by the foil within a cycle.

To quantify the aforementioned distance, we compare different length approximations of the path length ($\mathcal{L}$) covered by the TE in one non-dimensional period for a heaving foil. We estimate this length in three different ways:(a) step motion,(b) square wave and (c) sine wave. As shown in figure 6a (red curve) the step motion definition is equivalent to the use of $A$ to capture the covered length. As we see in figure 6b the square wave length $Sq$ captures $A$ in the vertical direction but also the streamwise length ($U_{\infty}/f$) traversed by TE in one period. Finally the trajectory length $Tr$ of the sine wave ( see figure 6c) captures the exact $\mathcal{L}$ traversed by the TE over the entire period:

[TABLE]

The utility of these three metrics is estimated via the agreement (collapse) of the neutral curves for different types of motions in the $f-\mathcal{L}$ domain. Here the $\mathcal{L}$ is normalised by $D$ and this is done so that $f$ can be expressed by $Sr$ at the same domain. Again figures 6a, 6b and 6c show the neutral lines for different kinematics in $\mathcal{L}_{D}-1/Sr$ domains for the aforementioned types of $\mathcal{L}_{D}$ namely: $A_{D},~{}\emph{Sq}_{D}$ and $\emph{Tr}_{D}$ respectively. The frequency is represented in the inverse form so that a power law between amplitude and frequency is depicted by a straight line. Since previous studies suggest that BvK reversal depends solely on $St_{A}$, we expect a collapse of the neutral line curves of different kinematics when plotted on an $A_{D}-1/Sr$ chart. However as seen in figure 6a these curves are heavily dependent on kinematics. On the other hand, as we move towards $Tr$ (which accurately represents the trajectory length) the collapse of the various neutral line curves is significantly improved (see figure 6c).

$Tr$ is suitable for characterising heave dominant kinematics since the motion of the entire chord can be represented just by the path traversed by the TE. This means that the displaced fluid per cycle of a heave dominant motion can be expressed as $Tr\cdot\mathcal{C}$. However, $Tr$ fails to capture the effects of a chordwise gradient $A$ present in pitch-dominated motions since most of the lateral motion is downstream of the pivot point for $\mathcal{P}<0.5$ (see figure 3b). This means that $Tr$ overestimates the displaced fluid for these particular cases and therefore it is still dependent on kinematics (see figure 6c). The effect of the pitching component can be incorporated into the metric by calculating the average trajectory length covered by the entire foil (chord) over one period:

[TABLE]

where, $s$ is the coordinate along the chord with $s=0$ at LE and $s=1$ at the TE.

The $D$-normalised chord averaged trajectory length $\mathcal{T}_{D}$ is plotted versus $1/Sr$ in figure 6d. The new metric demonstrates remarkable collapse of different kinematics on the curve corresponding to the neutral line of pure heave (where $Tr_{D}=\mathcal{T}_{D}$). Interestingly this curve follows the diagonal of a square and as we move towards more inviscid flows it can be expressed as $\mathcal{T}_{D}\cdot Sr\sim const=1$. More specifically for $Re=1173$ , $\mathcal{T}_{D}\cdot Sr\sim 1.035$ while for $Re=1173$ , $\mathcal{T}_{D}\cdot Sr\sim 1.015$. This product represents the average speed of the foil over one period. Thus the area in the lower right half of the straight line represents the zone where this average speed is less than $U_{\infty}$ while the the upper left is where the speed is larger than $U_{\infty}$. This provides a very simple physical interpretation of the $1/f-\mathcal{T}$ phase space where, in order for the foil to produce thrust, the kinematics has to be tuned in such a way that the chord averaged speed of the foil along the path over one period must be faster than the free stream. In other words, a path length based Strouhal number ($St_{\mathcal{T}}=f\mathcal{T}/U_{\infty}$) should be greater than 1. As the value of $St_{\mathcal{T}}$ increases further beyond 1, the wake of the foil becomes stronger producing higher and higher values of thrust.

The universality of the new length scale is evaluated by further examining the agreement of neutral line curves for various kinematic factors such as $f,~{}\mathcal{P},~{}\alpha$. Figure 5b shows all the kinematic options where a collapse was demonstrated. The agreement among these curves deteriorates for $Sr<0.12$ or $Sr>0.35$ perhaps due to the limitations of using two dimensional numerical simulations. This is consistent with the observations of Mittal & Balachandar (1995) who found that 2D simulations might result in large force fluctuations and erroneous wake patterns. Additionally non periodic wakes have been observed for pure pitch motion at $\mathcal{P}>0.6$ which could also contribute to the poor collapse. As $\alpha$ is essentially an indicator of the heave to the pitch ratio within a coupled motion case it has no real effect on periodicity at least up to $\alpha=20^{\circ}$ and thus the coupled motions presented here agree well with the pure pitch and heave test cases. Clearly the new metric sets a threshold for the drag-to-thrust wake transition of 2D flapping foils for the entire range of kinematics provided that the resultant wake is periodic.

3 Conclusions

Two dimensional simulations were conducted for a rigid flapping NACA0016 profile at low and high Reynolds numbers incompressible flow. The wake development towards a reversed BvK street was examined for a variety of harmonic motions, amplitudes and frequencies. At very low amplitudes and frequencies 2P wake patterns were observed for the coupled motions in agreement with either the pure heave or the pure pitch cases. This is due to the fact that at such low $Sr-A_{D}$ combinations the coupled motions are dominated by either the heaving or the pitching component.

In dimensionless Amplitude-Period maps various length scales were evaluated with respect to the neutral line of different motion types. It was revealed that the relationship between $A_{D}$ and the period is non-linear since the maximum distance from equilibrium cannot properly characterise the displaced volume (or area) required to overcome the drag forces. This is solved by calculating the length of the path traversed by the foil over one period of oscillation.

On a $1/Sr-Tr_{D}$ graph the neutral line of pure heave forms a linear curve $y(x)$ where $dy/dx=45^{\circ}$. Since $U_{foil}=Tr_{D}/(1/Sr)$ this means that thrust is achieved when $U_{foil}>U_{\infty}$. Furthermore the dimensionless chord average trajectories per cycle $\mathcal{T}_{D}$ of all motion types tested, collapse upon the pure heaving case. In other words the neutral lines of all test cases collapse on a trajectory-based Strouhal $St_{\mathcal{T}}\rightarrow 1$. Thus the new metric can serve as a universal length scale that captures the BvK reversal for every combination of harmonic two dimensional kinematics.

This novel method allows us to to parametrise drag-to-thrust wake transition of a simple two dimensional harmonically oscillating body via the chosen kinematics without the use of complex fluid dynamic equations. Knowing the onset of thrust for a flapping foil via a single parameter can significantly reduce the effort of designing sufficient biomimetic propulsors. Moreover it will enable scientists and engineers to describe and/or confirm observations regarding the thrust generating strategies and evolution of natural flyers and swimmers.

Acknowledgements

This research was supported financially by the Office of Naval Research award N62909-18-1-2091 and the Engineering and Physical Sciences Research Council doctoral training award [1789955].