Geometric Characteristics of Flapping Foils for Enhanced Propulsive Efficiency

Airfoils usually define cross-sectional profiles of wings, fins, and stabilizers of flying and swimming objects that are designed to create favorable lift-to-drag or thrust-to-drag ratios during locomotion. Prescribing a flapping motion to these foils is commonly seen in biological propulsors. Examples of this locomotion can be found in numerous natural species, including fish, birds, and insects that fly and swim [1]. A flapping motion allows for greater propulsive efficiency and maneuverability while creating less disturbance than traditional fixed or rotating foils [2]. Recently, studying flapping foils has supported the scientific community to pioneer new technologies in flow energy harvesters [3–5], micro-aerial vehicles [6,7], and bio-inspired underwater vehicles [8,9] in order to take advantage of various effective fluid–structure interaction-related mechanisms employed in nature. The recent focus on this field of research led to a large body of literature understanding how oscillating foils behave with varying kinematics, flow regimes, and a limited range of their geometric features [10,11]. However, how the geometric form of a flapping foil impacts flow physics, wake dynamics, and the resultant performance parameters has drawn limited attention from researchers. This work aims to provide a holistic view of foil shapes in a flapping motion, their effect on foil performance, and the underlying governing mechanisms that strongly influence propulsive performance.

To complete this study, it is of utmost importance to understand how variations in geometric configurations of foils influence the underlying flow physics and wake mechanics to produce optimum performance metrics. The governing mechanisms involve the formation and dynamics of complex vortices and their interaction with each other to change pressure variations around the foil, impacting its performance. In this context, Gopalkrishnan et al. [12] conducted a landmark experimental study for active vorticity control methods by introducing constructive and destructive interference between vortices produced around and over the surface of flapping foils. They examined force histories and vortex visualizations around flapping foils placed in the wake of a D-section cylinder at a Reynolds number of 550. Their work established correlations between performance metrics and complex vortex interference phenomena and confirmed the applicability of these results at a higher Reynolds number of 20,000 and a Strouhal number of 0.20. The study concluded that a destructive interference of the vortex shed from the cylinder with the trailing-edge vortex of the foil gave a higher propulsive efficiency. Contrarily, constructive interference with these coherent flow structures caused a higher thrust production.

Previous studies have been conducted investigating the foil thickness effect in flapping foils [13–16]. Ashraf et al. [13] numerically investigated the influence of varying thickness and camber on a foil, performing heaving and pitching simultaneously, for Reynolds numbers ranging from 200 to 200,000. They found that, in higher Reynolds number cases, a thicker leading-edge led to higher efficiencies and higher thrusts, and that camber had little to no impact on the overall performance of the foil. Other studies from Wang et al. [15] and Yu et al. [16] examined the effects of the thickness of a foil, the location of its maximum value along the chord length, and camber on the energy produced by a flapping foil energy extractor at a Reynolds number of 1200. In their study, they found that efficiency increased with increasing the maximum thickness up to 25% of the chord length. For a further increase in this parameter, the efficiency quickly dropped off. They highlighted that the overall performance could be significantly different under the same flow conditions and kinematics for different foil shapes. From these studies, we conclude that foil shape has a significant impact on the performance of the foil over a flapping motion. They also provide some preliminary insight into what shapes perform better than others. However, each of these studies rely on using a small number of standard NACA foil series to generate foils. Within each series, only camber and thickness are changed. Because of this, their studies are limited to observing how these parameters affect performance for that particular shape profile. To build upon this groundwork, this study aims to find the mechanisms behind changing shapes around the most optimal flapping foil shapes, while varying not only the thickness but the whole profile of the foil shape, allowing for a more universally applicable conclusion to be achieved.

A few studies have also varied from using standard foil series shapes in studying flapping foils. A flapping foil shape optimization study was completed by van Buren et al. [17] to create shapes of foils, which were then optimized experimentally and numerically. They investigated the performance of these foils at Reynolds numbers ranging from 1,000 to 10,000, with reduced frequencies of 0.4 and 1.0. They concluded that the foils exhibiting optimal performance appeared to have thicker leading edges. Recently, a study from Han et al. [18] numerically optimized the geometry of a flapping foil at Reynolds numbers, varying from 1000 to 10,000. They used multiple parameterization methods to generate different foil-like shapes, and reported that both a thicker leading-edge and a sharper trailing-edge were the primary contributors to improve the performance of the optimal foils. Finally, Kelly et al. [19] studied the effect of class-shape transformation (CST) parameters on the efficiency and thrust performance of a flapping foil, utilizing six CST coefficients in their study. They concluded that the second and third coefficients had the largest effect on the performance. These studies establish that a thicker leading edge and thinner trailing edge appear in flapping foils optimized for efficiency, breaking out from the standard foil series. However, they do not investigate the flow physics around changing the shape near this optimal shape regime.

While significant progress has been made in the study of flapping foils with a variety of cross-sectional shapes, a complete study on the hydrodynamics resulting from changing the shape through an extensive range of flapping foil shapes, breaking out from the traditional methods of using standard foil series shapes with varying only thickness, is required to fully understand the effect that shape has on performance in a flapping foil. In this work, we focus on conducting a holistic study of near optimal foil shapes undergoing a flapping motion. To focus the scope of the study, we look at foil shapes in the range of the optimal foils found in Ref. [18]. Particularly, we examine the influence the shape of a foil has on its propulsive performance and the governing physical mechanisms responsible for variations in performance parameters. For this purpose, we employ a CST parameterization method to generate geometric parameters that are varied to create unique shapes of foils. Then, we prescribe its flapping kinematics and flow conditions in the range of those observed in natural flying and swimming. Next, computational simulations are performed for flows over these newly generated foils undergoing oscillations using our in-house immersed boundary method-based direct numerical flow solver. From these simulations, we find the range of foil shapes that give optimal performance. By analyzing the resulting flow created by these simulations, we draw conclusions about the flow structures that lead to optimal performing foils and can be used to drive foil shape design in future applications.

In this section, we present the details of our kinematics, foil shape generation, immersed-boundary method-based computational solver, and associated boundary conditions.

where f is the frequency of the motion, $θ°$ is the maximum pitching angle, and, and $ϕ$ is the phase angle between the pitching and heaving motions.

Previously, van Buren et al. [20] analyzed the effect of changing $h°, θ°$⁠, and $ϕ$ on the performance of flapping foils. They concluded that the most efficient motion occurred at $ϕ=270 deg$⁠. In addition, they found that cases with $h°=0.375c$ produced the highest efficiencies. Furthermore, they found the best performance with respect to efficiency in cases with $θ°=15 deg$⁠. Therefore, these important result provides the basis for our choices of values for $ϕ, h°$⁠, and $θ°$ for our present work.

Next, the methodology for varying the shape of the foils is determined. Many previous studies focused on using standard shapes of foils, such as those belonging to NACA series. In order to ensure the capturing of a large spectrum of geometric configurations of foils, we need a parameterization method that allows us to fully control the shape of a foil. The use of such a method allows the complete exploration of the design space, rather than being limited by standard series foil shapes. Reviews of the previously developed and employed parameterization methods to create better performing foils were provided by various researchers [21,22]. In Ref. [21], multiple parameterization methods, including analytic method, singular value decomposition (SVD) method, class-shape transformation (CST) method, discrete method, domain element method, free form deformation (FFD) method, PARSEC method, PDE method, and Splines method, are compared. While this study did not provide a recommendation for a single method as superior, it does highlight SVD and CST for being more efficient in defining foil shapes with few parameters than the others. In another study [22], B-splines, CST, SVD, and PARSEC were analyzed in more details from the aspects of accuracy and efficiency. CST and SVD were also found to be the two most efficient methods in defining foil shapes with a less number of independent parameters.

where $xN1·(1−x)N2$ is the class function, $KjN·xj−1·(1−x)N−j$ is the shape function, $BjN$ are the basis functions, and a_j are the CST coefficients. First, in order to give the general foil shape of a round leading-edge and pointed trailing-edge, we define $N1=0.5$ and $N2=1.0$⁠. It is important to mention that the order of the Bernstein polynomial in the shape function is set to 5, which makes N = 6. This was determined from the work of Kulfan et al. [24] and earlier shown to be sufficient for 2D foil shapes [18,19,25]. This value of N implies that we have six parameters that are used next to determine the foil shape $a=(a1,a2,…,a6)$⁠, where a_j are the parameters corresponding to Eq. (3) and a is the vector used to denote all these parameters. Shown in Fig. 1(c) are the basis functions for the CST with $N1=0.5, N2=1.0$⁠, and N = 6. It can be observed that each function corresponds to different locations of peak values along the chord. This particular feature enables the manipulation of thickness in each region. It is made possible by increasing or decreasing the value of a coefficient with the basis function that peaks in the corresponding region. Figure 1(b) illustrates this process through an example, where each coefficient is increased about 10% from the NACA0012 value to vary the foil shape. We notice that using a larger coefficient values make the resulting foil thicker in the regions where its basis function is the largest. Additionally, the CST-generated NACA0012 foil, the values of which are given in Ref. [18], is shown alongside the actual NACA0012 foil. This confirms that the CST method accurately portrays this foil shape using N = 6.

where ν is the kinematic viscosity of the fluid. The solver discretizes Eq. (6) using a cell-centered collocated arrangement of the primitive variables and utilizes a Cartesian grid-based sharp interface immersed-boundary method. A second-order Adams-Bashforth scheme is employed to compute the convective terms, and the diffusion terms are computed using an implicit Crank–Nicolson scheme. This setup eliminates the need for remeshing as the foil undergoes oscillations. It was previously applied to simulations for flow around flapping foil successfully [19,26–28] and was also extensively used for a broad range of biological flows [29–35]. More details about the solver can be found in Ref. [36].

The Cartesian computational grid used for our current simulations is shown in Fig. 2(a). The domain size is $16c×14c$⁠, and the total grid contains about 954,000 (993 × 961) nodes. To accurately capture the foil geometries and flow boundary layers, an extremely dense region is assigned around the foil, with a minimum grid spacing of $Δmin=0.00227c$⁠. Outside this layer, a grid with resolution $Δmid=0.0125c$ is generated to resolve the wake features behind the flapping foil with sufficient accuracy. To exclude the effects of the grid on the hydrodynamic force calculation, grid independence study, shown in Fig. 2(b), is performed on three sets of grids with different minimum grid spacing, $ΔCoarse=0.00468c, ΔNominal=0.00227c$⁠, and $ΔFine=0.00175c$⁠, with a NACA0012 foil undergoing previously defined kinematics. From this grid study, the difference between the mean thrust coefficients obtained from the nominal and fine grids is 0.98%. Furthermore, the difference in the maximum thrust coefficients from the nominal and fine mesh is $1.15%$⁠. Hence, it is adequate to use nominal grid for our simulations.

To setup the flow in the computational domain, the left boundary is set to a velocity inlet with a constant incoming flow speed of $2.5c/T$⁠. A zero-gradient boundary condition is considered to the upper and lower sides, and a zero stream-wise gradient is prescribed to the outlet on the right side, as shown in Fig. 2(a). A homogeneous Neumann boundary condition is used for pressure at all boundaries.

where $CT¯$ and $CPw¯$ are the cycle-averaged coefficients of thrust and power, respectively.

where $A=2h0$ is the peak to peak trailing edge amplitude. These values are aligned with the previously studies [18,19] and lie within the operating range of natural swimmers [1,27].

In this section, we explain our findings for varying shapes of foils, their performance metrics, and the governing flow physics.

First, we perform simulations for flows over flapping NACA0012 foil generated with the CST technique shown in Fig. 1. To understand the impact of each coefficient on the foil's performance, we then change a single CST coefficient while holding the remaining ones for the conventional NACA0012 profile, resulting in a range of foils similar to Fig. 1. The shapes are varied to find the peak efficiency for each coefficient. In order to better compare the results to the NACA0012 foil the shapes originate from, the relative efficiency is shown. This is computed as $ηR=η¯ηNACA0012¯$⁠. Additionally, we utilize these coefficients from the optimized foil reported in Ref. [18] to ensure that the current range of coefficients include the ones with maximum η for each coefficient.

The results shown in Fig. 3(a) exhibit η for each case, normalized by the propulsive efficiency of a NACA0012 foil with the varying values of coefficients. It is important to highlight that a₅ and a₆ have a very narrow range of η, and are best around the values for a NACA0012 foil. This aligns with our hypothesis and previous studies showing that a geometric configuration with a thicker trailing-edge and without a thicker leading-edge is inefficient [18]. Looking at Fig. 1(b), increasing a₅ or a₆ without changing the other coefficients creates a thick trailing-edge. Similarly, the values of a₁ are optimal around those associated with a NACA0012 foils, but it offers a slightly broader range of values for high η. The underlying reason could be illustrated by considering Fig. 1(b), where an increasing a₁ too much makes the foil into a teardrop shape that is large only at its leading edge. In the plot, for a₄, we observe that there is a much larger range of values for higher η, and that the most efficient foil has a larger coefficient value than that for a NACA0012 foil. However, the overall effect of increasing this coefficient is small, and leads to only an enhancements of 2% in η at best. The foils corresponding to changes in a₃ are seen to have a significantly larger influence on the maximum η along with an increasing trend of η as a₃ increases gradually from that of a NACA0012 foil up to a steep drop-off point at an a₃ value of about 0.5. At the maximum, changing only this particular coefficient increases η by over 10% from that of a NACA0012 foil. Finally, variations in a₂ have the largest impact on η, because increasing its value from that of a NACA0012 foil improves this performance metric substantially, until a sharp drop is noticed for $a2=0.32$⁠. At its peak, changes in a₂ alone increases η by over 13%. These plots also reveal that the performance of a flapping foil is very sensitive to the geometric shape, meaning that very small variations in the shape of a foil lead to large differences in performance parameters. It is evident that η is particularly sensitive to thickness and shape of the foil between the locations $0.18c$ and $0.5c$ along the chord. It can be seen from Fig. 1, where the basis functions, corresponding to a₂ and a₃, are dominant in that region. Near the peaks of a₂ and a₃, η of the foil changes by more than 2% for every 1% change in a_i value. Furthermore, these results show that, for other coefficients, such as a₄ and a₅, the performance is not sensitive to variations in the foil's shape in that region.

These new findings challenge the assumption made by many previous studies related to flapping foils [37–40]. These research investigations focused on various aspects of flapping foils, including 3D effects, kinematics, and flow features etc., solely relying on standard shapes of foils and assuming the results not being sensitive to variations in their geometric configurations. Our present study reveals that small changes in the shape can give large changes in results, and that this effect varies for each coefficient. Therefore, simple variations in thickness of a flapping foil, belonging to a classical series of geometric configurations, do not necessarily scale to the whole domain of flapping foils.

Next, we examine a wider range of shapes of foils, changing multiple parameters sequentially. First, from this study, we consolidate our findings by observing the effect that the maximum thickness and maximum thickness location have on the performance of the foil. Contour plots are presented in Figs. 3(b) and 3(c) that provides a map of η and C_T as a function of δ_max and S_max. Selected cases with a constant S_max changing δ_max (⁠$Aδ, Bδ$⁠, and $Cδ$⁠) and constant δ_max changing S_max (A_S, B_S, and C_S) are chosen to present in depth analysis, each including the most efficient foil (⁠$Bδ,S$⁠). The efficiency values in this plot, shown by the color contour, are normalized by the NACA0012 foil's efficiency. The middle region of this contour colored in red exhibits a higher η. The most efficient propulsive performance are shown by the foils with $0.18<δmax<0.22$ and $0.26<Smax<0.33$⁠. A smaller width of the red region here compared to its height is also apparent, which hints that η is less sensitive to variations in S_max and more dependent on δ_max. Additionally, the average thrust coefficient (⁠$C¯T$⁠), normalized by the NACA0012 thrust coefficient, is plotted as a function of maximum thickness (δ_max) and maximum thickness location (S_max). From this plot, it can be seen that the best performing foils in thrust are generally centered in two different regions, one around maximum thickness of δ_max = 0.19 and max thickness location of S_max = 0.375, and the other around a maximum thickness of δ_max = 0.18 and max thickness location of S_max = 0.26. Unlike efficiency, the maximum thrust production occurs in two regions and is equally sensitive to the maximum thickness and the location of maximum thickness along the foil.

In order to isolate the effects on performance due to maximum thickness for further investigation, three foil shapes are chosen with a constant maximum thickness location to review in depth. The cases chosen are shown in Fig. 3 and details are given in Table 1. The S_max value of 0.3 was chosen to include the highest efficiency foil studied. In the table, it is observed that the most efficient foil, $Bδ$⁠, has a large increase in efficiency with a small increase in thrust compared to the thinner foil, $Aδ$⁠. The thicker foil, $Cδ$⁠, experiences a large drop in both efficiency and thrust.

To illustrate temporal variations in hydrodynamic performance parameters and investigate the performance differences further, temporal profiles of C_T and C_Pw over a single flapping cycle of motion are shown for the three aforementioned cases in Figs. 4(a) and 4(b). The coefficients of thrust and power are calculated using Eqs. (8) and (9). The figure also shown the shapes of the foils (c), along with the thrust produced (d) and power consumed (e) along the length of the foil over a cycle of motion. From the figure, we see that both the thrust and power exhibit their respective peaks at around $t/T=0.05$ and again at $t/T=0.55$⁠, which corresponds to the end of the heaving motion. C_T and C_Pw then dip back down, reaching their minima at $t/T=0.30$ and $t/T=0.80$⁠. This corresponds to the point immediately after the foil direction change occurs (⁠$t/T=0.30$⁠).

where p is the hydrodynamic pressure. In the vorticity, an obvious difference occurs in the V1 vortex on the top leading edge of the foil. In case $Aδ$⁠, the leading edge vortex begins to separate from the foil body, generating vortex V1 and resulting in a stronger low pressure area on the top front edge of the foil. Conversely, in case $Bδ$ the leading edge stays attached well to the body of the foil. There is however still a low pressure region formed on the top leading edge of the foil. In the thinner foil $Aδ$⁠, the normal to this low pressure region faces primarily in the $+y$ direction. Alternatively, in the thicker foil $Bδ$ the normal to this surface faces more in the –x direction, promoting the overall thrust. The larger force in $+y$ in $Aδ$ leads to the increased power consumption through this portion of the motion, and prevents the increased low pressure region from enhancing the thrust significantly over the value observed in foil $Bδ$⁠.

In comparing $Cδ$ to $Bδ$⁠, there is a significant decrease in thrust production that begins as the foil changes directions and exists through much of the next stroke (⁠$t/T=0.75$⁠). In comparing the countour plots of thrust production, a large low thrust region is observed in the back 60% of the foil is observed for foil $Cδ$ that does not appear for the other cases. In the vorticity and pressure plots at this time (Fig. 5), it is observed that case $Cδ$ has significantly more breakdown of the leading edge vortex along the back half of the foil, following the point of maximum thickness. Specifically, a significantly larger V2 vortex has built up and completely separated from the foil body. In the corresponding pressure contour, it can be seen that this also corresponds to a lack of a high pressure region between V1 and V2 observed for the other foils. This significant drop in pressure leads to the overall decrease in the net force in –x, giving the drop in cycle-average thrust observed. This separation occurs due to the larger thickness giving an overall steeper foil behind the maximum thickness point, promoting separation of the leading edge vortex. The vortex separation prevents constructive interference between the V2 vortex with the trailing edge vortex, as is seen in $Bδ$⁠. The lack of constructive interference leading to a lower thrust production is consistent with Ref. [12].

Overall, the thickness is found to be a sensitive parameter for efficiency due to the breakdown of the leading edge vortex, along with the proportion of the normal surface that points in the forward direction. In the case where the foil is thinner, there is less of a component in the thrust direction and an increase in horizontal force. While the thrust is compensated by the breakdown of the leading edge vortex at the front edge of thinner foils giving better suction force, this only increases the already higher horizontal force, consuming additional power. Because of this, increase in thickness gives the gradual increase in efficiency observed in Fig. 3(b). When the foil becomes too thick, however, the leading edge vorticies begin to break down behind the point of maximum thickness. This breakdown removes the high pressure region along the back of the foil that contributes significantly to thrust production. Because the transition from smoother vorticity to the total breakdown of the leading edge vortex is more abrupt, the efficiency drops quickly when the thickness becomes too high in Fig. 3. These results show that generally, a larger leading edge is advantageous, which is consistent with the results found in Ref. [13].

To similarly isolate the effects on performance due to maximum thickness location along the foil for further investigation, three foil shapes are chosen with a constant maximum thickness to review in depth. Detailed values for the cases chosen are given in Table 2. The δ_max value of 0.21 was chosen to include the highest efficiency foil studied. In the table, it is observed that the most efficient foil, B_S, has a large increase in efficiency with a small increase in thrust compared to the other two.

To investigate this further, the temporal profiles of C_T and C_Pw over a single flapping cycle of motion is shown for the three aforementioned cases in Figs. 6(a) and 6(b). The coefficients of thrust and power are calculated using Eqs. (8) and (9). The figure also shown the shapes of the foils (c), along with the thrust produced (d) and power consumed (e) along the length of the foil over a cycle of motion. The overall trends from the figure are similar to the maximum thickness cases, with both C_T and C_Pw exhibit their respective peaks at around $t/T=0.55$ and reaching their minima at $t/T=0.80$⁠.

In comparing case A_S with the maximum efficiency case B_S, a significantly reduced thrust produced is observed as the foil completes the change of direction, centered around $t/T=0.8$⁠. In the contour plot for thrust produced, a larger low thrust region is observed at this time around 50 percent of the foil chord length. To investigate the unsteady vortex dynamics further, the vorticity and pressure are plotted in Fig. 7. In this figure, a similar phenomena is observed to the previous section in the thickest foils. In case A_S, the leading edge vortex has broken down and separated completely from the body of the foil. Specifically, the two large vorticies on the top edge, V1 and V2, have completely separated from the body. The high pressure region along the foil body between V1 and V2 is greatly reduced compared to other cases, giving less thrust. The more abrupt transition from the beginning of the foil to the maximum thickness point in this case causes a higher velocity along the top edge behind the maximum thickness point, leading to the destabilization and breakdown of the leading edge vortex before it occurs for other cases.

In comparing case B_S with case C_S, a decrease in thrust produced and increase in power consumption is observed at the end of the heaving motion as the change in direction is initiated. In the contour plots for thrust and power, it is observed that the thrust difference occurs primarily at the front of the foil, around 25% of the chord, and the power difference occurs at the front of the foil body. In the pressure and vorticity plots shown in Fig. 7, a similar phenomena is observed to the thinnest foils discussed in the previous section. The leading edge vortex starts to break down near the front edge, leading to a larger vortex building up (V1). This gives a similar larger high pressure region causing an increase in the lateral force, resulting in an increase in power consumption. This can be attributed to the delay in the maximum thickness of the foil occurring so far back along the chord that it acts in the front as a thinner foil would. Additionally, a higher pressure region is observed developing behind the V1 vortex. This occurs still before the maximum thickness point, so the surface of the foil has a normal component in the –x direction. The higher pressure at this point gives the reduction in thrust produced for this region discussed above. This occurs because of the aforementioned geometric phenomena leading to the formation of V1, along with the maximum thickness point being far enough back along the foil to allow the higher pressure to be a detriment to the thrust production in the foil.

Similar to δ_max, efficiency is shown to be sensitive to S_max due to a breakdown of the leading edge vortex around the point of maximum thickness, which does not occur for the most efficient foil shape but happens for the others. The optimal values are seen around S_max = 0.3, which is consistent with the findings in Ref. [15].

Next, the effect of maximum thickness and maximum thickness location at different Reynolds numbers is studied. To accomplish this, select foil shapes are simulated at varying Reynolds numbers, to confirm the applicability of our results across varying flow regimes. First, cases varying the maximum thickness with a constant maximum thickness location along the foil are tested, changing the Reynolds number to 2000, 5000, and 10000. The resulting propulsive efficiencies are shown in Fig. 8(a). Next, the same process is repeated changing the maximum thickness location, keeping the maximum thickness constant. The results are shown in Fig. 8(b). First, an overall trend of higher Reynolds number giving higher efficiency is observed across all of the cases. In the cases changing maximum thickness, a similar curve is shown for all Reynolds numbers, though the peak efficiency values occur at smaller thicknesses for lower Reynolds numbers. These results are consistent with the findings of Ashraf et al. [13], showing that for as Reynolds number increases, the optimal thickness for efficiency increases also. A similar phenomena is observed to Fig. 5, however the complete leading edge vortex separation occurs in thinner foils at lower Reynolds numbers, giving the curve shape observed. In the maximum thickness location cases, a similar phenomena occurs with the same overall efficiency curves, however, the optimal thickness location for efficiency is shifted toward the back of the foil. Similar to the thickness cases, the leading edge vortex separation that makes the smallest maximum thickness location cases inefficiency occurs more in the lower Reynolds number cases. This gives the shift in efficiency plot when changing Reynolds number seen in the figure.

In this work, we conducted a holistic study of the shape of flapping foils. Specifically, we examined the effect of foil shape on performance and the underlying mechanisms that cause the variation in performance. We used a CST parameterization method to create geometric parameters that are varied to create unique foil shapes. We then prescribed a flapping motion and flow conditions in the range of the ones seen in biological swimming. Computational fluid dynamics simulations of these shapes in a flapping motion were completed using an in-house immersed boundary method-based direct numerical solver. From the initial study varying a single CST parameter at a time, we found that the efficiency of a foil was very sensitive to the shape of the foil. Specifically, we saw a significant change in performance when adjusting the thickness of the foil in the range of 18 to 50 percent of the chord-length. Outside of this range, the shape of the foil had little effect on the efficiency.

Next, we studied a larger subset of foil shapes for efficiency, and found that the maximum thickness and the location of the maximum thickness had an optimal range of 0.2c and 0.3c, , respectively. Sensitivity of efficiency to the maximum thickness was greater than it was to its location along the foil. Upon analyzing the continuous coefficients of thrust and power, along with the vortex structures and pressure fields for varying each parameter individually, we found that better attachment of the leading edge vortex, leading to constructive interference with the trailing edge vortex. Additionally, a favorable pressure gradient during the turning from the up-stroke to the down-stroke was observed, which lead to the increase in thrust generated and decrease in power consumed, , respectively, for the most efficient cases. Similar results were shown in the optimal of both the maximum thickness and maximum thickness location, with variations in either leading to poorer adhesion of the leading edge vortex, either in the front of the foil before the point of maximum thickness in cases with thinner foils (V1), or at the back of the foil behind the maximum thickness point in cases with thicker foils (V2).

Overall, we found that a thicker leading-edge on a foil led to more favorable pressure gradients to lower power consumption during a flapping motion, while a thicker foil results in a drop in thrust produced due to the lack of adhesion of the leading edge vortex preventing constructive interference with the trailing edge vortex. We saw that not only is the overall performance of a flapping foil highly sensitive to its shape, but also the vortex structures behind the flapping performance is significantly different depending on the choice of the foil shape. Because of this sensitivity and change based on the shape of a foil, care should be taken when selecting foil shapes to use for bio-inspired underwater robotics, and careful shape selection can lead to significant performance benefits. Our results are contrary to the assumption made in many previous studies that performance and vortex dynamics results in flapping foils will scale similarly when changing the shape of the foil. It may not be sufficient to use a standard series foil shapes when alternate foil shapes have different performance and flow features.

This work was supported by ONR MURI N0014-14-1-0533, NSF CNS-1931929, and the NSF NRT program. JK would also like to acknowledge ASME FED GSS for their support in this work. Additionally, thanks to the University of Virginia Research Computing Group for the availability of the Rivanna supercomputing cluster.

• National Science Foundation (Award ID: CNS-1931929; Funder ID: 10.13039/100000001).

• Office of Naval Research Global (Award ID; Funder ID: 10.13039/100007297).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.