3.1 Results from full simulations

The results presented in this section were obtained from all 400 tandem-foil simulations, and the thrust force was calculated directly by taking the integral of the pressure and shear stress on the surface of the foil, although there is minimal contribution from the latter. To allow comparison to the virtual hind foil results presented below in § 3.2, the current results are referred to as ‘measured’ and represented by the additional subscript ‘ $m$ ’.

Thrust. Figure 2 shows the normalised thrust coefficient contour plots for the hind foil $C_{T,h,m}^{\ast }$ and it can be seen that the hind foil of the tandem-foil arrangement is able to produce up to twice the thrust of the single foil under certain parametric combinations. This shows that the hind foil of the tandem system can be much more effective than a single foil, and would potentially make very effective propulsion systems. However, certain parametric combinations lead to lower thrust, so the geometry and kinematics of the system must be accurately controlled to obtain the benefit. As shown by the diagonal bands of high and low thrust, the thrust force is highly dependent on both phase and spacing. These have the form of parallel ‘ridges’ and ‘valleys’. This has previously been observed for single frequencies (Boschitsch et al. Reference Boschitsch, Dewey and Smits2014; Gong et al. Reference Gong, Jia and Xi2015, Reference Gong, Jia and Xi2016), but until now this effect has never been observed over a range of frequencies. The fact that such a large thrust augmentation occurs for all these Strouhal numbers shows that the effect is important over a wide range of conditions.

Figure 2.

Contour plots of the measured mean thrust coefficient of the hind foil with respect to the single foil ( $C_{T,h,m}^{\ast }$ ), for all phase ( $\unicode[STIX]{x1D711}$ ) and spacing ( $S$ ) combinations over four Strouhal numbers, each calculated using 80 tandem-foil simulations. Red is high thrust and blue is low thrust, white represents the value of a single foil.

Efficiency. The contours of propulsive efficiency of the hind foil (figure 3) have the same form of diagonal bands as the thrust coefficient contours (figure 2). For the phase–spacing combinations that lead to higher thrust the efficiencies are generally the same as the efficiency of a single foil. Where the hind foil produces a lower thrust, the efficiencies are much lower than a single foil.

Figure 3.

Contour plots of the measured mean efficiency of the hind foil with respect to the single foil ( $\unicode[STIX]{x1D702}_{h,m}^{\ast }$ ), for all phase ( $\unicode[STIX]{x1D711}$ ) and spacing ( $S$ ) combinations over four Strouhal numbers, each calculated using 80 tandem-foil simulations. Red is high efficiency and blue is low efficiency, white represents the value of a single foil.

Qualitative flow description. To determine the flow characteristics that give rise to the high- and low-performance contour bands, two representative cases were chosen for further analysis. For a Strouhal number of 0.4 and a spacing of $S=2$ , figure 2 shows that the highest performance occurs at ( $\unicode[STIX]{x1D711}=7\unicode[STIX]{x03C0}/4$ ) and the lowest performance occurs at ( $\unicode[STIX]{x1D711}=3\unicode[STIX]{x03C0}/4$ ). These two cases will hereafter be referred to as the ‘high-performance’ and ‘low-performance’ cases respectively.

Figure 4.

Contour plots of instantaneous vorticity magnitude ( $\unicode[STIX]{x1D6FE}$ ) and time-averaged streamwise velocity field ( $\overline{u}$ ) for the single, high-performance $(S=2,\unicode[STIX]{x1D711}=7\unicode[STIX]{x03C0}/4)$ and low-performance case $(S=2,\unicode[STIX]{x1D711}=3\unicode[STIX]{x03C0}/4)$ , all at a Strouhal number of 0.4. For the vorticity fields the fore foil is midway through the upstroke, clockwise vorticity is indicated with blue and anticlockwise with red. In all figures, the vertical grid lines show the locations of the quarter-chord points of the fore and hind foils, and the horizontal grid lines show the midline of the domain and the upper and lower limits of heave motion. Videos of the vorticity fields for the single, high-performance and low-performance cases can be found in the on-line supplementary movies and are entitled movie 1, movie 2 and movie 3 respectively.

Figures 4(a), 4(c) and 4(e) show the instantaneous vorticity fields for the single, high-performance and low-performance cases respectively. These are snapshots of movies that can be found on-line in the supplementary movies as movie 1, movie 2 and movie 3 available at https://doi.org/10.1017/jfm.2017.457. Figures 4(b), 4(d) and 4(f) are the time-averaged streamwise velocity fields for the single, high-performance and low-performance cases respectively. All of these are for a Strouhal number of 0.4. The typical flow structure of a single thrust-producing flapping foil is known to consist of two primary vortices of alternating sign shed by the foil per cycle (Gopalkrishnan et al. Reference Gopalkrishnan, Triantafyllou, Triantafyllou and Barrett1994; Anderson et al. Reference Anderson, Streitlien, Barrett and Triantafyllou1998; Read et al. Reference Read, Hover and Triantafyllou2003). This wake structure is similar to a drag-producing von Kármán street of a bluff body, although the signs of the vortices on each side of the ‘street’ is reversed in the flapping foil case because it produces thrust instead of drag. The instantaneous vorticity field of the single foil case (figure 4 a) shows that these primary vortices exist along with smaller secondary vortices which appear as the shear layers roll up. When averaged over time this arrangement of vortices induces a thrust-producing jet behind the foil (figure 4 b).

Figure 5.

Velocity vectors for flapping foil with different inlet conditions, showing approximate location of vortex, and induced velocity components.

The vorticity field for the high-performance tandem case (figure 4 c, where $S=2$ and $\unicode[STIX]{x1D711}=7\unicode[STIX]{x03C0}/4$ ) shows that the hind foil weaves between the vortices that are shed from the fore foil. This increases the vorticity on the front surface of the hind foil, and therefore also the strength of the vortices that are shed into the wake. These strong vortices become interspersed with the vortices from the fore foil, pairing up with them and creating a ‘double von Kármán street’. This creates a faster velocity behind the foils than the single foil case (figure 4 d), which accounts for the higher thrust. This type of vortex–foil interaction that produces high forces on the hind foil has been referred to as a ‘high-thrust’, ‘constructive’ or ‘coherent wake’ mode in the literature.

The vorticity field of the low-performance tandem case (figure 4 e, where $S=2$ and $\unicode[STIX]{x1D711}=3\unicode[STIX]{x03C0}/4$ ) shows that the hind foil encounters (or comes very close to) the vortices shed by the fore foil. This decreases the strength of the vorticity on the surface of the hind foil, and therefore the strength of the vortices shed in the wake. A more dispersed and weaker jet is present behind the foils in this case than the high-performance tandem case. This type of interaction has previously been called the ‘low-thrust’, ‘destructive’ or ‘branched wake’ mode.

In summary, if the hind foil weaves between the vortex street of the fore foil then its thrust will be increased. Conversely, if the hind foil encounters these vortices then its thrust will be attenuated. Inspection of a variety of Strouhal numbers, phases and spacings (not shown here) confirmed that for each frequency the high- and low-performance cases have the same type of vortex interaction, regardless of the particular phase and spacing at which these cases occur. Two more important observations are, first, that the high- and low-performance cases are always separated by a phase difference of $\unicode[STIX]{x03C0}$ , and second, that the time-averaged jet velocity behind the fore foil is generally unaffected by the phase and spacing, implying that the vortex advection speed between the foils is constant for each Strouhal number.

The effects of the vortex-induced velocities ( $u_{\unicode[STIX]{x1D6FE}}$ and $v_{\unicode[STIX]{x1D6FE}}$ ) on the lift and thrust production in figure 5 shows why the hind foil must have a phase and spacing which cause it to weave between the vortices of the fore foil for the highest thrust production. It can be reasoned that if the crossflow component of the vortex-induced velocity ( $v_{\unicode[STIX]{x1D6FE}}$ ) is in the same sense as the heave velocity of the foil ( ${\dot{h}}(t)$ ) then the magnitude of the induced velocity ( $U_{I,h}$ ), and therefore the thrust production at that instant in time will be increased. Also, if $(v_{\unicode[STIX]{x1D6FE}}/u_{\unicode[STIX]{x1D6FE}})>({\dot{h}}(t)/U_{\infty })$ then the angle of attack ( $\unicode[STIX]{x1D6FC}_{h}$ ) will increase.

To determine when these conditions occur for the hind foil, the instantaneous streamwise and crossflow velocity fields are considered (figures 6(a) and 6(b) respectively). The streamwise velocity behind the fore foil is always higher than the freestream so will always act to increase the induced velocity of the hind foil regardless of its specific kinematics. This is in agreement with the jet-type wake in figure 4(b). On the other hand, the crossflow velocity alternates between every primary vortex, so to increase its induced velocity, the hind foil would need to move in opposite direction to this crossflow velocity, which is why the hind foil must weave between the incoming vortices for the highest thrust production. These observations indicate that the thrust and efficiency of the hind foil can be calculated using the velocity field in the wake of the fore foil.

Figure 6.

Instantaneous vortex-induced velocity fields for the single foil at $St=0.4$ .

3.2 Virtual foil method

To move beyond the qualitative analysis of the previous section, we performed a quantitative quasi-steady analysis of a ‘virtual’ hind foil to determine the reasons for the performance variations of the hind foil. Quasi-steady analyses have been used to predict the forces on flapping foils to a reasonable degree of accuracy for a range of situations (Jensen Reference Jensen1956; Weis-Fogh Reference Weis-Fogh1972, Reference Weis-Fogh1973; Ellington Reference Ellington1984a ,Reference Ellington b ; Sane & Dickinson Reference Sane and Dickinson2002; Nakata, Liu & Bomphrey Reference Nakata, Liu and Bomphrey2015), so is a suitable method of analysis for the current problem.

The performance variation observed above must be due to the hind foil encountering different incoming flow conditions due to the wake of the fore foil at each phase–spacing combination. Therefore we postulate that a quasi-steady-state analysis of the properties of the inflow condition of the hind foil will elucidate the fundamental parameters that control the performance augmentation and also yield a reasonable approximation of the thrust and efficiency of the hind foil of figures 2 and 3. The virtual foil method developed predicts the thrust and efficiency of a hind foil by using the instantaneous velocity field behind a single flapping foil, which must be known at many instances throughout the flapping cycle. This velocity field is used along with the equations of motion of a virtual hind foil to determine the induced velocity magnitude and angle of attack of this foil. Simple steady-state aerodynamics are then used to calculate the instantaneous lift, drag and moment at each time, from which are calculated the thrust, side force, power and efficiency of the virtual hind foil. An analysis of the accuracy of this virtual foil method of predicting the performance of a single foil is given in appendix C. It shows that the virtual foil method analysis can predict the cycle-averaged thrust and efficiency of a single foil with reasonable accuracy and is therefore suitable for the tandem flapping foil system. The virtual foil method has four main steps:

Step 1: Measure inflow conditions of virtual hind foil.

The velocity components due to the fore foil are determined from the instantaneous velocity field of a simulation of a single foil (these could also be determined from experiments). At each time the position of the virtual hind foil is calculated from its kinematics (2.1). The inflow velocities ( $u_{\unicode[STIX]{x1D6FE}}$ and $v_{\unicode[STIX]{x1D6FE}}$ ) at this foil position are found from the velocity field obtained from the simulations of a single foil.

Step 2: Calculate induced velocity and angle of attack of virtual hind foil.

Induced velocity. For the case of a single or fore foil, which encounters a steady incoming flow, the induced velocity ( $U_{I,f}$ ) consists of the vectorial addition of the freestream velocity, $U_{\infty }$ , and the heave velocity of the foil, ${\dot{h}}(t)$ , as shown in figure 5. For the case of the hind foil of a tandem system which encounters an unsteady incoming flow including the vortex street of the fore foil, the induced velocity ( $U_{I,h}$ ) consists of the freestream and the heave velocities as before, but now with the additional velocity components caused by the vortices ( $u_{\unicode[STIX]{x1D6FE}}$ and $v_{\unicode[STIX]{x1D6FE}}$ ) (figure 5).

Angle of attack. For the single or fore foil, the angle of attack can be calculated using the heave velocity ( ${\dot{h}}(t)$ ), the freestream velocity ( $U_{\infty }$ ) and the pitch angle of the foil ( $\unicode[STIX]{x1D703}$ ). For a heave to pitch phase of $\unicode[STIX]{x1D6F9}=\unicode[STIX]{x03C0}/2$ , the angle of attack of the fore foil is

For a foil in an unsteady flow with the additional velocity components caused by the vortices ( $u_{\unicode[STIX]{x1D6FE}}$ and $v_{\unicode[STIX]{x1D6FE}}$ ) this becomes

The normalised square of the induced velocity ( $U_{I,h}^{2\ast }$ ) and angle of attack ( $\unicode[STIX]{x1D6FC}_{h}^{\ast }$ ) of the virtual foil show dependence on phase and spacing (figure 7 a,b), displaying similar contours to the measured thrust augmentation of figure 2. This confirms that both of these quantities affect the mean thrust production of the hind foil. The variation of the angle of attack is much larger than the variation of the induced velocity, which shows that the dominant parameter affecting the thrust augmentation of the hind foil is the angle of attack. This is observed for the entire range of Strouhal numbers considered, although only the data for $St=0.4$ is shown here for brevity. This novel analysis shows that the thrust augmentation of in-line tandem flapping foils is primarily due to an alteration of the angle of attack of the hind foil to the modified inflow velocity angle in the presence of the incoming vortex street.

Figure 7.

Mean inlet conditions of the hind foil calculated using the virtual foil method for all phase–spacing combinations at a Strouhal number of 0.4, and normalised by the corresponding values of a single foil. (a) The normalised mean of the square of the induced velocity $(U_{I,h}^{2\ast })$ and (b) the normalised mean of the angle of attack $(\unicode[STIX]{x1D6FC}_{h}^{\ast })$ .

Step 3: Calculate the lift, drag and moment of the foil using steady-state aerodynamic theory.

Now that the instantaneous induced velocities and angles of attack have been calculated, they can be used to determine the lift, drag and moment at each instance during the flapping cycle. Ignoring any dynamic effects, the instantaneous lift and drag forces ( $L$ and $D$ ) and moment ( $M$ ) on a steady foil can be found using

where $\unicode[STIX]{x1D70C}$ is the density of the fluid. The coefficient of lift is directly proportional to the angle of attack of the foil ( $\unicode[STIX]{x1D6FC}$ ), and can be approximated by

The coefficient of drag and for the NACA 0016 profile was found using a panel method to be

The moment coefficient was assumed to be linear up to $\unicode[STIX]{x1D6FC}\approx 10^{\circ }$

Step 4: Calculate the thrust and efficiency.

Figure 8.

Contour plots of the mean thrust coefficient of the hind foil with respect to the single foil ( $C_{T,h,v}^{\ast }$ ), for all phase ( $\unicode[STIX]{x1D711}$ ) and spacing ( $S$ ) combinations for four Strouhal numbers, each calculated using 1 single foil simulation using the virtual foil method. Red is high thrust and blue is low thrust, white represents the value of a single foil.

The lift, drag and moment determined in the previous step can be used to calculate the thrust coefficient and efficiency of the virtual hind foil. The drag and lift forces act parallel and perpendicular to the induced velocity angle, which is ( $\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D703}$ ). The thrust and side forces, which are parallel and perpendicular to the freestream velocity are therefore

The thrust coefficient, power and efficiency can then calculated using (2.3a – c ), (2.5) and (2.4) respectively.

Figures 8 and 9 show the thrust coefficient and efficiency contours for four Strouhal numbers calculated using the virtual foil method outlined above. These results were obtained from one single foil simulation for each frequency, and to differentiate them from the measured results, they are referred to as ‘virtual’ and denoted with the additional subscript ‘ $v$ ’. The results given here are normalised by their respective values of the thrust of the single foil calculated using the virtual foil method rather than by the single foil value from the measured results. This normalisation avoids offset errors due to the overestimation of the mean thrust coefficients of the single foil as seen in figure 14(a) in appendix C. These virtual foil contours are very similar to the measured thrust coefficient contours of figure 2, which shows that the virtual foil method is reasonable at reproducing the trends observed in the full measured results. However, there are some differences, which are discussed in the following section.

Figure 9.

Contour plots of the efficiency of the hind foil with respect to the single foil ( $\unicode[STIX]{x1D702}_{h,v}^{\ast }$ ), for all phase ( $\unicode[STIX]{x1D711}$ ) and spacing ( $S$ ) combinations for four Strouhal numbers, each calculated using 1 single foil simulation using the virtual foil method. Red is high efficiency and blue is low efficiency, white represents the value of a single foil.

3.3 Comparison of virtual foil results to measured results

The differences between the mean thrust coefficient contours for the measured and virtual foil results must be quantified to determine the limitations of the virtual foil method. This has been achieved by characterising the thrust contours using four simple parameters. Consider figure 10, which is a representation of the contours of figures 2 and 8. As previously discussed, the thrust augmentation bands are parallel ridges (high thrust) and valleys (low or no thrust). If we consider a high-thrust ridge, a straight line can be fitted to the locus of maxima (at every spacing) along the ridge using least-squares. Only two parameters are required to describe this line: the slope and intercept. Additionally, the average value of the high thrust along the ridge and low thrust along the valley characterises the amount of thrust augmentation achieved along those bands at each Strouhal number. These four characteristic parameters for the measured contours can be compared to the values obtained from the virtual foil method to determine its ability at predicting the thrust of the hind foil.

Figure 10.

Schematic of the characterisation of the thrust augmentation bands using the slope and intercept of the high-thrust contour ridge.

The intercept of the high-thrust contour with the phase axis is the ‘zero-spacing phase’ and is represented by $\unicode[STIX]{x1D711}_{0}$ . It can be interpreted as the phase shift that is required to maximise the thrust coefficient of the hind foil when the spacing is zero. For the measured results, $\unicode[STIX]{x1D711}_{0}$ decreases exponentially with Strouhal number, and tends to zero with increasing frequency as seen in figure 11(a). This implies that at a spacing of zero and high frequency, the phase that maximises the thrust coefficient of the hind foil is zero. In this situation the foils are effectively moving together and acting as one foil. The virtual foil results on the other hand do not show this trend, and the intercept remains close to zero for all frequencies. It is likely that this is because the virtual foil method calculates the thrust of the foil instantaneously, and does not account for any time delay or advection effects of the flow.

Figure 11.

Characteristic values of contour plots for measured and virtual foil data: (a) intercept $\unicode[STIX]{x1D711}_{0}$ against St, (b) intercept $\unicode[STIX]{x1D711}_{0}$ against $\sqrt{C_{T}}$ , (c) slope $\text{d}\unicode[STIX]{x1D711}/\text{d}S$ against St, (d) minima $\overline{C_{T,min}^{\ast }}$ and maxima $\overline{C_{T,max}^{\ast }}$ against St. Measured data (calculated using figure 2 and represented by circles) and virtual foil method data (calculated from figure 8 and represented by triangles). Lines are solid for measured data and dashed for virtual foil method.

A more meaningful metric with which to compare the intercept is the characteristic velocity in the wake of the fore foil. Intuitively, the phase offset that is observed can be interpreted as the time taken for the vortices shed by the fore foil to ‘interact’ with the hind foil including the time taken by the vortices to reach the hind foil and advect past it. Given that the vortices shed by the fore foil have a range of sizes and strengths, it might not be feasible in every case to explicitly measure the time taken for this interaction. Therefore, a proxy for this advection speed is necessary. The square root of the thrust coefficient ( $\sqrt{C_{T}}$ ) is used as it represents the characteristic velocity of the jet in the wake of a thrust-producing foil. Figure 11(b) shows that $\unicode[STIX]{x1D711}_{0}$ is linear with $\sqrt{C_{T}}$ , as $\unicode[STIX]{x1D711}_{0}$ reduces with increasing value of the thrust of the fore foil. Therefore, as the velocity of the jet in the wake of the fore foil increases (with increasing thrust), the advection speed of the vortices increases and hence the time taken for these vortices to interact with the hind foil reduces. Although the virtual foil method is unable to predict the intercept of the contour bands accurately, the linear relationship of the measured results can be used to determine the necessary phase shift with which to adjust the results from the virtual foil method.

The slope of the thrust coefficient contour bands ( $\text{d}\unicode[STIX]{x1D711}/\text{d}S$ ) is the change in phase ( $\text{d}\unicode[STIX]{x1D711}$ ) over a given change in spacing ( $\text{d}S$ ) (figure 10) and is a measure of the relative importance of the phase with respect to spacing. The value of $\text{d}S$ is the primary component which determines the slope of the contours, as it represents the characteristic advective wake speed of the vortices generated by the fore foil. Figure 11(c) plots the slope of the contour bands for all Strouhal numbers for both the measured and virtual foil results. There is good agreement between the measured data and the virtual foil results, showing that the virtual foil method is able to predict the slope of the contour bands well. The slope of the contours becomes steeper with increasing Strouhal number, and this trend is linear over the current Strouhal number range. This means that at higher Strouhal numbers, the wake speed is faster, so the hind foil interacts with more vortices in a certain time interval. This is consistent with the interpretation of the phase offset discussed previously.

The mean value of the high thrust contour is calculated by taking the average of the maximum thrust coefficient for each spacing, and the minimum is determined in a corresponding fashion. The minimum thrust decreases with Strouhal number, showing that the detrimental foil–vortex interactions are stronger at higher frequencies. This trend is captured well with the virtual foil method. However, the virtual foil method overestimates the maximum thrust for most frequencies (figure 11 d), which can be explained by examining the efficiency of the single foil. Figure 14(b) shows that the efficiency of the single foil increases sharply up to the maximum at $St=0.3$ , after which it decreases more slowly. This is the inverse trend of the maximum thrust coefficient for the measured results of figure 11(d) and shows that when the single foil has a high efficiency, the maximum thrust augmentation of the hind foil is lower. We propose that this is because when the foil is flapping with kinematics that give it a high efficiency, it is already operating close to its maximum efficiency, so further increases in thrust are proportionally harder to obtain.

Knowledge of the slopes, intercepts, maxima and minima of the phase–spacing contour plots over a range of frequencies is powerful for the design of tandem flapping foil propulsors as it allows prediction of the parametric combinations that lead to high thrust without measurements of the entire parameter space. The contour plot of any Strouhal number between the frequencies given here can be reconstructed by using the slope and offset values given in figures 11(c) and 11(b).

For the three lowest Strouhal numbers (figure 9 ( $St=0.25$ not shown)), the virtual foil method predicts the form of the efficiency contours, but overpredicts the efficiency over the whole phase–spacing range. For certain phase–spacing combinations for the highest two Strouhal numbers, the virtual foil method calculates a power close to zero, and therefore predicts efficiency values that unrealistically low or high. We consider these errors to be because the virtual foil method does not account for additional losses in the flow due to the hind foil, and more work is required to understand this fully.

The virtual foil method not only leads to insights regarding the mechanism of the thrust augmentation of tandem flapping foils, and about the upstream inflow condition of immersed bodies, but it also greatly reduces the number of experiments required to investigate the time-averaged thrust over the entire spacing–phase–frequency parameter space of tandem flapping foils. This method shows that instead of running a tandem case of every combination of these parameters, just one run of a single foil is all that is required to predict the performance augmentation.