1. Introduction
The interplay between fins shapes, elastic properties, hydrodynamic forces and passive and controlled kinematics is the subject of persistent and active research (Videler Reference Videler1975; Geerlink & Videler Reference Geerlink and Videler1986; Tangorra et al. Reference Tangorra, Davidson, Hunter, Madden, Lauder, Dong, Bozkurttas and Mittal2007; Lauder Reference Lauder2015; Puri et al. Reference Puri, Aegerter-Wilmsen, Jaźwińska and Aegerter2018). Replicating the flexible fins or the body of fish using hydrofoils with controlled motion programs has proven a powerful tool to investigate the kinematics and propulsive forces of swimmers (Shelton, Thornycroft & Lauder Reference Shelton, Thornycroft and Lauder2014). Methods to quantify the fluid velocity field are typically based on imaging tracer particles seeded into the flow corresponding to particle imaging velocimetry (PIV) or particle tracking velocimetry (PTV) (Maas, Gruen & Papantoniou Reference Maas, Gruen and Papantoniou1993; Dracos Reference Dracos1996; Raffel et al. Reference Raffel, Willert, Wereley and Kompenhans1998; Pereira et al. Reference Pereira, Stüer, Graff and Gharib2006). The vortex wakes and thrust production of flapping foils with various geometries and flexibilities have been extensively characterized in the literature (Triantafyllou, Techet & Hover Reference Triantafyllou, Techet and Hover2004; Godoy-Diana, Aider & Wesfreid Reference Godoy-Diana, Aider and Wesfreid2008; Bohl & Koochesfahani Reference Bohl and Koochesfahani2009; Kim & Gharib Reference Kim and Gharib2010; Green, Rowley & Smits Reference Green, Rowley and Smits2011; David et al. Reference David, Jardin, Braud and Farcy2012; Marais et al. Reference Marais, Thiria, Wesfreid and Godoy-Diana2012; Shinde & Arakeri Reference Shinde and Arakeri2014; David, Govardhan & Arakeri Reference David, Govardhan and Arakeri2017; Lucas, Dabiri & Lauder Reference Lucas, Dabiri and Lauder2017; Muir, Arredondo-Galeana & Viola Reference Muir, Arredondo-Galeana and Viola2017). Furthermore, many researchers have resorted to particle velocimetry experiments in order to quantify the flow field of aquatic animal appendages and bioinspired synthetic fins (Blickhan et al. Reference Blickhan, Krick, Zehren, Nachtigall and Breithaupt1992; Stamhuis & Videler Reference Stamhuis and Videler1995; Müller et al. Reference Müller, Van Den Heuvel, Stamhuis and Videler1997; Drucker & Lauder Reference Drucker and Lauder1999; Lauder Reference Lauder2000; Müller, Stamhuis & Videler Reference Müller, Stamhuis and Videler2000; Drucker & Lauder Reference Drucker and Lauder2001; Müller et al. Reference Müller, Smit, Stamhuis and Videler2001; Nauen & Lauder Reference Nauen and Lauder2002a,Reference Nauen and Lauderb; Drucker & Lauder Reference Drucker and Lauder2005; Müller & Van Leeuwen Reference Müller and Van Leeuwen2006; Tytell Reference Tytell2006; Lauder & Madden Reference Lauder and Madden2007; Tangorra et al. Reference Tangorra, Davidson, Hunter, Madden, Lauder, Dong, Bozkurttas and Mittal2007; Müller, van den Boogaart & van Leeuwen Reference Müller, van den Boogaart and van Leeuwen2008; Tytell, Standen & Lauder Reference Tytell, Standen and Lauder2008; Tangorra et al. Reference Tangorra, Lauder, Hunter, Mittal, Madden and Bozkurttas2010; Flammang et al. Reference Flammang, Lauder, Troolin and Strand2011a,Reference Flammang, Lauder, Troolin and Strandb; Dewey, Carriou & Smits Reference Dewey, Carriou and Smits2012; Esposito et al. Reference Esposito, Tangorra, Flammang and Lauder2012; Ren et al. Reference Ren, Hu, Wang and Wen2016a,Reference Ren, Yang, Wang and Wenb; Mwaffo et al. Reference Mwaffo, Zhang, Romero Cruz and Porfiri2017).
The research field making use of PIV/PTV-based pressure evaluation covers a large collection of hydrodynamics problems, from micro-channel junction flow and turbine blades to animal locomotion (Gresho & Sani Reference Gresho and Sani1987; Baur & Köngeter Reference Baur and Köngeter1999; Gurka etal. Reference Gurka, Liberzon, Hefetz, Rubinstein and Shavit1999; Jakobsen, Dewhirst & Greated Reference Jakobsen, Dewhirst and Greated1999; Fujisawa et al. Reference Fujisawa, Nakamura, Matsuura and Sato2006; Liu & Katz Reference Liu and Katz2006; van Oudheusden, Scarano & Casimiri Reference van Oudheusden, Scarano and Casimiri2006; Murai et al. Reference Murai, Nakada, Suzuki and Yamamoto2007; van Oudheusden et al. Reference van Oudheusden, Scarano, Roosenboom, Casimiri and Souverein2007; van Oudheusden Reference van Oudheusden2008; Windsor Reference Windsor2008; Jardin, David & Farcy Reference Jardin, David and Farcy2009; Lorenzoni et al. Reference Lorenzoni, Tuinstra, Moore and Scarano2009; Khodarahmi et al. Reference Khodarahmi, Shakeri, Sharp and Amini2010; de Kat & van Oudheusden Reference de Kat and van Oudheusden2012; Ragni, van Oudheusden & Scarano Reference Ragni, van Oudheusden and Scarano2012; de Kat & Ganapathisubramani Reference de Kat and Ganapathisubramani2013; van Oudheusden Reference van Oudheusden2013; Panciroli & Porfiri Reference Panciroli and Porfiri2013; Dabiri et al. Reference Dabiri, Bose, Gemmell, Colin and Costello2014; Joshi, Liu & Katz Reference Joshi, Liu and Katz2014; Tronchin, David & Farcy Reference Tronchin, David and Farcy2015; Lucas et al. Reference Lucas, Dabiri and Lauder2017; McClure & Yarusevych Reference McClure and Yarusevych2017; Mwaffo et al. Reference Mwaffo, Zhang, Romero Cruz and Porfiri2017). Accurate and non-invasive methods to measure the fluid forces directly on the surface of flapping fins are essential to investigate the mechanisms of underwater propulsion. In the present work, we combine fluid velocity measurements from three-dimensional (3-D) particle tracking velocimetry with hydrodynamic stress calculations based on the Navier–Stokes equation to obtain well-resolved dynamic maps of fluid forces on the surface of flexible fin-like foils operating in different flow regimes. We address the long-standing question of how the geometry and kinematics of oscillating hydrofoils influence their spatio-temporal distributions of hydrodynamic forces. We focus on the effects of flapping frequency, width, length and edge curvature (straight versus bilobed). This approach offers a framework to investigate how specific morphological and kinematic features can constitute an advantage for certain hydrodynamic functions. The implications of this study reach into the fields of animal aquatic locomotion, as well as the engineering design of biomimetic underwater vehicles. Such propulsive systems are gaining importance due to their potential applications in the monitoring, maintenance and exploration of underwater environments (Low & Willy Reference Low and Willy2006; Lauder et al. Reference Lauder, Anderson, Tangorra and Madden2007; Low Reference Low2007; Zhou & Low Reference Zhou and Low2012; Qing Ping et al. Reference Qing Ping, Wang, Dong, Shang and Tan2013).
2. Material and methods
2.1. Particle tracking velocimetry
We use a three-dimensional, three-component particle tracking velocimetry approach to quantify the flow generated by the synthetic fins. Technical aspects of three-dimensional PTV are well described in Pereira et al. (Reference Pereira, Stüer, Graff and Gharib2006). The basic principle of PTV relies on tracer particles seeded in the flow, illuminated by a laser beam and imaged at regular time intervals to track the position of each particle, allowing the subsequent reconstruction of the fluid velocity field. Because the aim is to calculate the hydrodynamic stress tensor from the velocity fields, this experimental method presents several benefits compared to other techniques such as tomographic or scanning stereo-PIV, the most outstanding one being the possibility to instantaneously capture the whole flow volume, allowing for the straightforward reconstruction of the 3-D velocity vectors everywhere inside that volume.
The different components of the experimental set-up are shown in figure 1 with their technical specifications. The working fluid is water and the seeding tracers are polyimide particles with a diameter of 
$\sim$50 
$\mathrm {\mu }$m, which have been used in previous particle velocimetry experiments to investigate the flow fields generated by fish fins (Flammang etal. Reference Flammang, Lauder, Troolin and Strand2011a). The measurement volume in our experiments is 
$50\ \textrm {mm}\times 50\ \textrm {mm}\times 20\ \textrm {mm}$, with an approximate number of 
$8\times 10$
$^{4}$ seeding particles inside that interrogation domain. Three cameras mounted on a plate in a triangular arrangement are used to record the 3-D fluid velocity vectors 
$\boldsymbol {u}=(u_{x},u_{y},u_{z})$. The flow chamber has transparent walls on three sides and a fixation wall on one side for inserting the fin model, actuated with a servomotor fixed outside the chamber. The hydrofoils are mounted horizontally inside the tunnel, with their rotation axis parallel to the 
$z$ axis and perpendicular to the cameras plate. A frontal perspective would make the detection of the particles more arduous in front of the fin surface. Water tanks are connected to both sides of the flow chamber in a recirculating system, and a pump pushes water inside a pipe from one tank to the other to control the upstream fluid velocity. A flow straightener is installed at the inlet of the tunnel to insure laminar incoming stream. All the experiments were performed using the V3V-9800 system (TSI Incorporated), which is characterized in Lai et al. (Reference Lai, Pan, Menon, Troolin, Castaño-Graff, Gharib and Pereira2008). The laser double pulse timing and the cameras capturing frequency are synchronized in a method called frame straddling: the time difference between a pair of position fields is determined by the time difference between a pair of laser pulses (
$\delta t = 2.5$ ms), and the rate at which velocity fields are recorded (80 Hz) is half the camera acquisition frequency, yielding a time separation of 
${\rm \Delta} t=12.5$ ms between the velocity fields. The reconstruction of the particles positions is based on a 2-D Gaussian fit of the particles intensity distributions and a triplet search algorithm is used for the 3-D positions fields. The velocity vectors are computed by tracking the particle displacements between subsequent laser pulses, using a relaxation method to achieve a probability-based matching. These processing steps are conducted using the V3V software (version 2.0.2.7). An appropriate combination of median filter, velocity range and smoothing filter was applied to the raw velocity fields to reduce the noise level and remove the worst outliers. Moreover, a mask is applied over the hydrofoil during image processing to avoid detecting ghost particles inside its boundary. Finally, the raw vectors are interpolated on a regular 3-D grid using Gaussian interpolation, yielding a final spatial resolution of 0.75 mm in each direction. Temporal and spatial resolutions were chosen along recommendations offered in previous studies where the flow velocity field was used to compute hydrodynamic pressure (de Kat & van Oudheusden Reference de Kat and van Oudheusden2012; Dabiri et al. Reference Dabiri, Bose, Gemmell, Colin and Costello2014; Wang et al. Reference Wang, Gao, Wei and Wang2017). The grid points inside the hydrofoil boundary, where no particle is detected owing to the mask, are not attributed a velocity value and are not involved in the hydrodynamic stress calculation. The closest grid points to the real fin boundary define a virtual boundary, where the surface distributions of hydrodynamic forces are evaluated.
Experimental set-up: (1)camera triplet (4 MP, 85 mm lenses, sensor size 
$11.3\ \textrm {mm}\times 11.3\ \textrm {mm}$, magnification 0.3, max. frequency 180 Hz) arranged in a triangular configuration on a plate located at 
$\sim$465 mm from the middle of the (2) flow chamber, illuminated by a (3) laser beam (double-pulsed, Nd:YAG, wavelength 532 nm, max. 120 mJ pulse
$^{-1}$) expanded by a pair of cylindrical lenses and deflected by (4) a mirror. (5) Recirculating system composed of two tanks (total volume of 
$\sim$160 l) and a pump connected to a pipe. (6) Sketch of the flow chamber with inner dimensions (in mm), top view (
$x$–
$z$ plane) and frontal view (
$x$–
$y$ plane), measurement volume (
$\sim 50\ \textrm {mm}\times 50\ \textrm {mm}\times 20\ \textrm {mm}$) indicated with red rectangles, enclosing the trapezoidal fin. (7) Example of an image captured by one of the three cameras, with tracer particles and hydrofoil midline clearly visible (direction of the free-stream indicated by blue dashed arrows and direction of the normal and tangential forces indicated by green arrows).
The foil midline is visually tracked over time in the 
$x$–
$y$ plane (see label 7 of figure 1). For each time frame, points are manually superimposed on the fin midline and fitted with a polynomial of degree 2. The 3-D hydrofoil is reconstructed based on that fitted quadratic curve, assuming that deflection occurs only along one axis. The virtual object has a larger width than the actual foil (0.725 mm away from the real surface on each side) owing to the fact that the particles cannot be resolved directly at the fluid–solid interface. In order to evaluate if the reconstructed foil boundary is sufficiently close to the real surface where we want to extract the hydrodynamic pressure, we used the criterion that the virtual surface should be located within the fluid boundary layer. Based on the assumption of a laminar boundary layer, its thickness can be expressed as (in the simplified case of a flat plate, with U being the fluid velocity far from the plate and 
$\nu$ being the kinematic viscosity of the fluid) (Prandtl Reference Prandtl1952)
\begin{equation} \delta=5 \sqrt{\frac{\nu x}{U}}. \end{equation} We could estimate that the thickness of the fluid boundary layer is 
$>$1 mm everywhere on the fin except at the most proximal region (
$<$10 % of the fin length). The virtual surface thus remains inside the boundary layer of the real object. Even though the boundary layer of an oscillating fin involves more complex phenomena such as turbulence, separation and reattachment (Obremski & Fejer Reference Obremski and Fejer1967; Kobashi & Hayakawa Reference Kobashi and Hayakawa1980; Arnal Reference Arnal1984; Incropera et al. Reference Incropera, Dewitt, Bergman and Lavine2007; Kunze & Brücker Reference Kunze and Brücker2011), these flow effects would tend to make the boundary layer extend further away from the solid surface. Therefore, the assumption of laminarity of the boundary layer gives a conservative estimate of the boundary layer width and thus allows us to verify that the virtual foil boundary provides an accurate representation of the real surface distributions. This point was verified in a previous study (Dagenais & Aegerter Reference Dagenais and Aegerter2019), where we used a control volume analysis to compute the forces generated by a flapping fin in similar conditions, and showed that the results compared favourably to the integrated force distributions over the fin surfaces.
2.2. Morphologies and kinematics of the fin models
The synthetic fins are illustrated in figure 2 and characterized in table 1 and figure 3. Taking shape A
$_{1}$ as the reference geometry, A
$_{2}$ constitutes a shorter version and shape B, a wider version. Shape C has the same length as shape A
$_{2}$, but presents a bilobed trailing edge. The aspect ratio (AR) of a fin is defined as the square of the span (width at largest point) divided by the area. Even though the fin models do not mimic any fish species in particular, their geometries can be compared to caudal fins with low aspect ratios (
$\textrm {AR}=0.84$ for shape A
$_{1}$, 1.03 for shape A
$_{2}$, 1.67 for shape B and 1.1 for shape C). Examples include the Schistura genus, the Oryzias genus, the zebrafish (Danio rerio) and the platy fish (Xiphophorus maculatus), as shown in figure 3 (Sambilay Reference Sambilay2005; Offen et al. Reference Offen, Blum, Meyer and Begemann2008; Parichy et al. Reference Parichy, Elizondo, Mills, Gordon and Engeszer2009; Naruse, Tanaka & Takeda Reference Naruse, Tanaka and Takeda2011; Plongsesthee, Beamish & Page Reference Plongsesthee, Beamish and Page2012; Bohlen et al. Reference Bohlen, Petrtyl, Petra and Chhouk2016; Kottelat Reference Kottelat2017).
Four fin geometries under study. Principal axes of the fin are shown in red: proximo-distal axis (from 0 to 
$L$) and dorso-ventral axis (from 
$V$ to 
$D$).
Black dots: aspect ratio and rigid Strouhal number (
$St_r$, calculated based on the rigid projection of the peduncle) for each fin model. Horizontal lines: caudal fins aspect ratios for various fish species, derived from Sambilay (Reference Sambilay2005), Offen et al. (Reference Offen, Blum, Meyer and Begemann2008), Parichy et al. (Reference Parichy, Elizondo, Mills, Gordon and Engeszer2009), Naruse et al. (Reference Naruse, Tanaka and Takeda2011), Plongsesthee et al. (Reference Plongsesthee, Beamish and Page2012), Bohlen et al. (Reference Bohlen, Petrtyl, Petra and Chhouk2016) and Kottelat (Reference Kottelat2017). Dashed blue zone: range of Strouhal numbers typically associated with efficient propulsion for flapping foils (Triantafyllou et al. Reference Triantafyllou, Triantafyllou and Yue2000; Taylor et al. Reference Taylor, Nudds and Thomas2003). Dashed pink zone: range of Strouhal numbers for the adult zebrafish caudal fin, derived from Parichy et al. (Reference Parichy, Elizondo, Mills, Gordon and Engeszer2009), Palstra et al. (Reference Palstra, Tudorache, Rovira, Brittijn, Burgerhout, van den Thillart, Spaink and Planas2010) and Mwaffo et al. (Reference Mwaffo, Zhang, Romero Cruz and Porfiri2017).
Parameters of the fin models: length (
$L$), width at the tip (
$d$, distance between lobes tips for shape C), and Reynolds number (
$Re$).
To produce the foils, a rigid cast of their negative form was 3-D printed, then liquid PDMS (polydimethylsiloxane) was poured inside the cast and the supporting rod was inserted at the base. The material was cured for 36 h at 58 
$^{\circ }$C. The resulting flexible membranes have a thickness of 0.55 mm. The properties of the cured PDMS can be found in the MIT material property database (http://www.mit.edu/6.777/matprops/pdms.htm). Most importantly, the mass density matches that of water. A cantilever deflection set-up was employed by S. Puri (University of Zurich) to characterize the elastic properties (see Puri Reference Puri2018; Puri et al. Reference Puri, Aegerter-Wilmsen, Jaźwińska and Aegerter2018 for details), yielding a value of 0.8 MPa for the Young's modulus. An external velocity of 
$u_{\infty }=55$ mm s
$^{-1}$ is imposed in all experiments, and the fins are actuated at their leading edge with a sinusoidal pitching motion with angular amplitude 
$\theta _{0}=11^{\circ }$ in all cases. Shape A
$_{1}$ was selected to test the effects of pitching frequency, using 
$f_{1}=1.9$ Hz, 
$f_{2}=2.8$ Hz and 
$f_{3}=3.7$ Hz. The natural frequency of this fin model (first mode) was estimated by releasing it from a rest position and measuring the decaying oscillation of its tip (in water). A value of 2.4 Hz was found, in between the lowest and intermediate frequencies tested. For the experiments involving shapes A
$_{2}$, B and C, frequency 
$f_{2}=2.8$ Hz is used.
Two dimensionless parameters allow us to characterize the flow regime of the hydrofoils. The Strouhal number (
$St$) encompasses the propulsion dependence on tail oscillation. The Reynolds number (
$Re$) describes the viscous versus inertial effects and determines the transition from a laminar to a turbulent flow. These parameters are based on the kinematic viscosity of water (
$\nu$), the external fluid velocity (
$U=u_{\infty }$), the foil length (
$L$), the flapping frequency (
$\,f$) and the tip excursion amplitude (
$A$)
\begin{equation} Re=\frac{U L}{\nu},\quad St=\frac{f A}{U}. \end{equation} Caudal fins found in nature are extremely diverse and cover a vast range of flow parameters, depending on species but also on developmental stage and behaviour. The swimming velocities catalogued in the literature usually correspond to the maximal speeds which fish can sustain for a short time only, yielding 
$Re$ from 
$10^{4}$ to 
$5\times 10^{6}$ (using the animal body length) (Wardle Reference Wardle1975; Yates Reference Yates1983; Vogel Reference Vogel1994). These high values are not representative of slower gaits at which fish can also operate in natural conditions (Bainbridge Reference Bainbridge1960; Weerden, Reid & Hemelrijk Reference Weerden, Reid and Hemelrijk2013). For example, based on the vast literature about zebrafish hydrodynamics, it can be estimated that this species operates at 
$Re$ in the range [390–2600] (using the caudal fin length to allow direct comparison with isolated hydrofoils) (Parichy et al. Reference Parichy, Elizondo, Mills, Gordon and Engeszer2009; Palstra et al. Reference Palstra, Tudorache, Rovira, Brittijn, Burgerhout, van den Thillart, Spaink and Planas2010; Mwaffo et al. Reference Mwaffo, Zhang, Romero Cruz and Porfiri2017). In the present study, 
$Re$ is situated within [1062–1375] (table 1). Although not characteristic of fast bursts or maximal velocity locomotion, this flow regime can be compared to the case of fish cruising at approximately 0.5 to 0.9 body lengths per second, pertinent for certain species (Sambilay Reference Sambilay2005) and specific behaviours such as foraging, chemotaxis and exploration.
Furthermore, studies have demonstrated that fundamental flow features of flapping fin propulsion can be captured by experiments and simulations performed at 
$Re$ lower than values typically measured for fast swimming fish, down to of the order of 
$10^{3}$ (Lauder etal. Reference Lauder, Madden, Hunter, Tangorra, Davidson, Proctor, Mittal, Dong and Bozkurttas2005; Bozkurttas et al. Reference Bozkurttas, Dong, Mittal, Madden and Lauder2006; Buchholz & Smits Reference Buchholz and Smits2006; Kern & Koumoutsakos Reference Kern and Koumoutsakos2006; Mittal et al. Reference Mittal, Dong, Bozkurttas, Lauder and Madden2006; Bozkurttas et al. Reference Bozkurttas, Mittal, Dong, Lauder and Madden2009; Liu et al. Reference Liu, Ren, Dong, Akanyeti, Liao and Lauder2017). Indeed, for sufficiently high 
$Re$ (
$\geq 10^{3}$), this parameter plays a minor role compared to 
$St$ in defining the flow topology (Lentink Reference Lentink2008; Green et al. Reference Green, Rowley and Smits2011).
In the case of flexible foils, 
$St$ depends on the tip excursion amplitude and has to be measured rather than fixed prior to the experiments. Therefore, we additionally define the rigid Strouhal number (
$St_r$), an input parameter based on the excursion amplitude of an equivalent rigid fin (same length and angular amplitude at the base). For the six experimental cases, we obtain 
$St_r$ in the range [0.33–0.64] (figure 3). Strouhal numbers between 0.2 and 0.4 usually yield the highest propulsion efficiency (Triantafyllou, Triantafyllou & Yue Reference Triantafyllou, Triantafyllou and Yue2000; Taylor, Nudds & Thomas Reference Taylor, Nudds and Thomas2003). Nevertheless, wider ranges of 
$St$ within [0.2–0.7] were reported in various fish species and for different developmental stages (Eloy Reference Eloy2012; Eloy Reference Eloy2013; Weerden et al. Reference Weerden, Reid and Hemelrijk2013; Xiong & Lauder Reference Xiong and Lauder2014; Van Leeuwen, Voesenek & Müller Reference Van Leeuwen, Voesenek and Müller2015; Link et al. Reference Link, Sanhueza, Arriagada, Brevis, Laborde, González, Wilkes and Habit2017). For example, adult zebrafish flap their caudal fins with 
$St$ in a window of [0.37–0.52], depending on their swimming mode (Parichy et al. Reference Parichy, Elizondo, Mills, Gordon and Engeszer2009; Palstra et al. Reference Palstra, Tudorache, Rovira, Brittijn, Burgerhout, van den Thillart, Spaink and Planas2010; Mwaffo et al. Reference Mwaffo, Zhang, Romero Cruz and Porfiri2017).
Aside from their pertinence in the fish world, the experimental parameters (AR, 
$Re$ and 
$St$) are similar to those selected in previous research about fluid–structure interactions of flapping foils (Lauder et al. Reference Lauder, Madden, Hunter, Tangorra, Davidson, Proctor, Mittal, Dong and Bozkurttas2005; Godoy-Diana et al. Reference Godoy-Diana, Aider and Wesfreid2008; Dai et al. Reference Dai, Luo, Ferreira de Sousa and Doyle2012; Marais et al. Reference Marais, Thiria, Wesfreid and Godoy-Diana2012; Dewey et al. Reference Dewey, Boschitsch, Moored, Stone and Smits2013; Quinn, Lauder & Smits Reference Quinn, Lauder and Smits2014; Shinde & Arakeri Reference Shinde and Arakeri2014; Quinn, Lauder & Smits Reference Quinn, Lauder and Smits2015; David et al. Reference David, Govardhan and Arakeri2017; Floryan et al. Reference Floryan, Van Buren, Rowley and Smits2017; Liu et al. Reference Liu, Ren, Dong, Akanyeti, Liao and Lauder2017; Rosic et al. Reference Rosic, Thornycroft, Feilich, Lucas and Lauder2017; Zhu et al. Reference Zhu, Wang, Lewis, Zhu, Dong, Bart-Smith, Wainwright and Lauder2017). We explore a specific portion of the flow parameter space to understand how propulsive efficiency and surface distributions of hydrodynamic forces depend on shape and frequency at relatively low Reynolds numbers, a topic of high interest in the field of fin propulsion and for the practical design of underwater vehicles relying on bio-inspired undulating membranes.
2.3. Calculations
Based on the PTV-velocity fields, the total hydrodynamic stress tensor 
$\boldsymbol {s}$ is calculated at each node in the 3-D domain, which is the sum of the scalar pressure field 
$p$ (multiplied with the identity matrix and a factor 
$-$1) and the viscous stress tensor 
$\boldsymbol {\tau }$
\begin{equation} s_{ij}={-}p\delta_{ij}+\tau_{ij}. \end{equation} The viscous stress tensor depends only on the spatial derivatives of the velocity field (where 
$\mu$ is the dynamic viscosity)
\begin{equation} \tau_{ij}=\mu\left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right). \end{equation} Although it was included in the calculation, the viscous stress is much smaller than the pressure in the present case. Indeed, the fins operate in the inertial flow regime (
$Re>1000$), although theoretically very close to the transitional range 
$300<Re<1000$ (McHenry & Lauder Reference McHenry and Lauder2005), where the normal stress component (dominated by the pressure) is typically larger than the viscous tangential stresses by at least two orders of magnitude.
The pressure evaluation is based on the Navier–Stokes equation (Whitaker Reference Whitaker1968; Aris Reference Aris1990)
\begin{equation} \rho \frac{\textrm{D}\boldsymbol{u}}{\textrm{D}t}={-}\boldsymbol{\nabla}p+\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau}+\rho \boldsymbol{g}. \end{equation} The last term on the right corresponds to any type of body force such as gravity; it is included in the pressure term and omitted in the remaining development. The left side of the equation contains the material acceleration 
$D/Dt$. This term is evaluated in the Lagrangian frame, namely in the reference frame of the advected particle (Dabiri et al. Reference Dabiri, Bose, Gemmell, Colin and Costello2014). Each component of the pressure gradient can be calculated from the Navier–Stokes equation
\begin{equation} \frac{\partial p}{\partial x_i}= \left(-\rho \left( \frac{\partial}{\partial t}+\sum_j{u_j \frac{\partial}{\partial x_j}}\right) +\mu \sum_j{\frac{\partial^2}{\partial x_j^2}} \right) u_i. \end{equation} Instantaneous pressure fields are reconstructed through direct integration of the above equation with appropriate boundary conditions. We performed the pressure calculation using the queen2 algorithm from Dabiri et al. (Reference Dabiri, Bose, Gemmell, Colin and Costello2014), available at http://dabirilab.com/software/. A null pressure value is assumed on the external boundary of the domain, in the undisturbed flow. The validity of that assumption relies on the conditions 
$H/D \geqslant 2$, where 
$H$ is half of the domain size (Dabiri et al. Reference Dabiri, Bose, Gemmell, Colin and Costello2014). Substituting the tip excursion amplitudes for the characteristic dimension 
$D$, we conclude that our experiments lie just above that limit. For each node inside the domain, the pressure gradient is integrated along eight different paths (horizontal, vertical or diagonal) originating on the outside contour. The median from the eight resulting pressure values is finally selected. This algorithm offers the advantage of reasonable computation time even for large domains. The crucial step in the pressure calculation lies in the determination of the material Lagrangian acceleration. In the queen2 algorithm, a so-called pseudo-tracking scheme is applied. The position of a particle at an instant 
$t_{i+1}$ is approximated based on its initial location 
$\boldsymbol {x}_p(t_{i})$ and the velocity evaluated at this location, averaged between instants 
$t_{i}$ and 
$t_{i+1}$. The velocity of the particle at its estimated forward position at time 
$t_{i+1}$ is then employed to evaluate its acceleration at time 
$t_{i}$. The pressure gradient is evaluated in a quasi-3-D manner, using in-plane velocity derivatives only, but combining integration paths in both the 
$x$–
$y$ and the 
$x$–
$z$ planes. The details of the calculation are presented in Dabiri et al. (Reference Dabiri, Bose, Gemmell, Colin and Costello2014).
The total hydrodynamic stress tensor (2.3) is projected on the surface of the solid object to obtain a stress vector, which is the total force per unit area generated by the fluid on the foil. Each 
$i$-component of the hydrodynamic stress vector acting on a surface with a unit normal vector 
$\hat {e}_{n}=(n_x,n_y,n_z)$ (oriented outwards) is expressed as
\begin{equation} S_i(\boldsymbol{x},t,\hat{n})=s_{ij}(\boldsymbol{x},t) \cdot n_j. \end{equation} In § 3, two types of distributions are presented: (i) instantaneous forces (per unit length), obtained by the integration of the surface stress maps along the dorso-ventral axis, at five selected time points covering half a period of oscillation, and (ii) period-averaged stress maps, where the normal stress is averaged either in absolute value or with its sign. In each experiment, the flow field is collected over three flapping periods and 30 pairs of velocity fields (10/period) are selected for the evaluation of hydrodynamic forces. Instantaneous distributions are averaged over six equivalent time frames (considering the symmetry between the left and right strokes). The period-averaged distributions are based on the 30 time frames. The instantaneous force (per unit length) is decomposed into the 
$x$ and 
$y$ directions, yielding the thrust 
$F_x$ and lateral force 
$F_y$ acting on the fin. The thrust force corresponds to the useful power spent by the fin (when it points in the negative 
$x$ direction, propelling the system upstream). The lateral forces correspond to the wasted power. Hence, we define the efficiency ratio as the power invested by the fin into useful thrust divided by the total rate of work done on the fluid
\begin{equation} \eta = \frac{|{\langle F_x \rangle}| \langle \overline{u}_{\textit{fin}} \rangle}{|{\langle F_x \rangle}| \langle \overline{u}_{\textit{fin}} \rangle + \langle |{F_y}| \rangle 2Af}. \end{equation} This definition is equivalent to the Froude efficiency: it measures the ability of the fin to convert work into upstream propulsion (Eloy Reference Eloy2013; Quinn et al. Reference Quinn, Lauder and Smits2015). The over lines indicate a spatial average and the brackets indicate a time average over a period. Note that 
$F_x$ is averaged with its sign; the absolute value is applied after averaging, whereas 
$F_y$ is averaged in absolute value directly because the fin is wasting energy due to lateral forces from both sides. The fluid velocity is averaged over the foil virtual boundary rather than over the whole volume, in order to better capture the specific ability of each hydrofoil to generate downstream fluid motion (see table 2). This classification based on 
$\eta$ is compared to the usual efficiency classification based on the Strouhal number, where the window of efficient propulsion is considered to be [0.2,0.4] (Triantafyllou et al. Reference Triantafyllou, Triantafyllou and Yue2000; Taylor et al. Reference Taylor, Nudds and Thomas2003).
Kinematics and flow regime of the hydrofoils: phase lag between peduncle and tip angles (
${\rm \Delta} \phi$, in fraction of a period), total amplitude covered by the fin tip over a period (
$A$), streamwise velocity component averaged over three periods (
$\overline {u}_{\mathit {vol}}$, averaged over the whole volume, and 
$\overline {u}_{\mathit {fin}}$, averaged over the virtual surface enclosing the fin), Strouhal number (
$St$) and efficiency ratio (
$\eta$).
Many authors have calculated the denominator of (2.8) for pitching fins (input power) as the period-averaged product of torque and angular velocity at the base (Dewey et al. Reference Dewey, Boschitsch, Moored, Stone and Smits2013; Lucas et al. Reference Lucas, Thornycroft, Gemmell, Colin, Costello and Lauder2015, Reference Maas, Gruen and Papantoniou2017; Quinn et al. Reference Quinn, Lauder and Smits2015; Egan, Brownell & Murray Reference Egan, Brownell and Murray2016; David et al. Reference David, Govardhan and Arakeri2017; Floryan et al. Reference Floryan, Van Buren, Rowley and Smits2017; Rosic et al. Reference Rosic, Thornycroft, Feilich, Lucas and Lauder2017; Zhu et al. Reference Zhu, Wang, Lewis, Zhu, Dong, Bart-Smith, Wainwright and Lauder2017). This efficiency metric is appropriate for hydrodynamic experiments where forces and torque sensors are placed at the attachment rod of the flapping propulsor. In the present study, no load cells are involved and all information about the forces imparted by the foil on the fluid are extracted solely from the flow velocity field. In this context, a definition of 
$\eta$ based on local force components integrated on the fin surfaces is more suitable. A similar calculation using surface integration of hydrodynamic stresses was used in Liu etal. (Reference Liu, Ren, Dong, Akanyeti, Liao and Lauder2017).
The uncertainties on the pressure and force distributions can be obtained using error propagation arguments and an analysis of the inaccuracy on the Lagrangian path reconstruction (used in the material acceleration evaluation). The uncertainties on the particles positions in our PTV experimental set-up are 
$\sigma _x=\sigma _y \simeq 3.6\ \mathrm {\mu } \textrm {m}$ and 
$\sigma _z \simeq 32\ \mathrm {\mu }\textrm {m}$, which imply velocity uncertainties of 
$\sigma _{u_x}=\sigma _{u_y} \simeq 0.002$ m s
$^{-1}$ and 
$\sigma _{u_z} \simeq 0.018$ m s
$^{-1}$. Noise propagation from the velocity field to the material acceleration and to the integrated pressure field has been the object of many studies (Liu & Katz Reference Liu and Katz2006; Violato, Moore & Scarano Reference Violato, Moore and Scarano2011; de Kat & Ganapathisubramani Reference de Kat and Ganapathisubramani2013; van Oudheusden Reference van Oudheusden2013; Wang et al. Reference Wang, Gao, Wei and Wang2017). We can derive an expression for the pressure uncertainty, relevant for the present calculation scheme. As a first step, we express the material acceleration uncertainty (with 
$a_i=\textrm {D}u_i/\textrm {D}t$) as a sum of the error propagated from the velocity field (first term under the square root) and the uncertainty from the Lagrangian path line reconstruction (second term)
\begin{equation} \sigma_{a_i} = \sqrt{\frac{2\sigma_{u_i}^2}{{\rm \Delta} t^2} + \left( (\boldsymbol{\sigma}_u \boldsymbol{\cdot} \boldsymbol{\nabla})u_i \right)^2}. \end{equation} The pressure field is the result of a spatial integration, its uncertainty therefore depends on the spatial resolution 
${\rm \Delta} x_i$ and on the number of nodes 
$n$ crossed along the integration path. Moreover, the pressure integration algorithm involves a median polling among a collection of 
$N=8$ paths, which reduces the uncertainty further by a factor of 
$\sqrt {{{\rm \pi} }/{2(N-1)}}$ (Kenney & Keeping Reference Kenney and Keeping1962, pp. 32–35, 52–54, 211–212). Because any of the 
$x$, 
$y$ or 
$z$ directions can be followed by the integration paths, the estimated pressure uncertainty is based on the average of errors in all three directions (Dagenais & Aegerter Reference Dagenais and Aegerter2019)
\begin{equation} \sigma_p = \frac{\rho}{3} \sqrt{\frac{\rm \pi}{14}}\sum_{i=1}^{3}n {\rm \Delta} x_i \sqrt{\frac{2\sigma_{u_i}^2}{{\rm \Delta} t ^2} +\left( (\boldsymbol{\sigma}_u \boldsymbol{\cdot} \boldsymbol{\nabla})u_i \right)^2}. \end{equation} The local pressure uncertainties are propagated to the force uncertainties using classic error propagation through an integration step, and the resulting values are shown with error bars in § 3. Due to temporal and spatial averaging, these uncertainties are reduced by an additional factor of 
$\sqrt {N_{\mathit {time}} \times N_{\mathit {spatial}}}$; 
$N_{\mathit {time}}$ is the number of time frames employed in the period average (6 or 30 for the instantaneous and period averaged distributions, respectively); 
$N_{\mathit {spatial}}$ is 1 for the instantaneous distributions and 20 in the case of the period-averaged curves of figure 11, which are spatially averaged over the left/right and dorsal/ventral symmetric sides of the fin, as well as over five dorso-ventral rows for each curve.
3. Results
3.1. Hydrofoil kinematics
The geometry and flapping frequency directly affect the fluid–structure interactions of the hydrofoils. The deflection profiles of the midlines along the course of a flapping cycle are shown in figure 4. The midline excursions are used to evaluate 
$A$, the tip amplitude, which is needed to calculate 
$St$. The phase lag between the peduncle and the tip (
${\rm \Delta} \phi$) is determined by measuring the angle between the fin midline and the horizontal plane (at the peduncle and at the tip) at regular time points during the motion cycle. The results are listed in table 2. The largest phase lag is found for the wide fin geometry (shape B), with 
${\rm \Delta} \phi$ twice as large as for shape A
$_{1}$. This tip recoil is associated with a smaller excursion amplitude at the tip (
$A$). Increasing the flapping frequency of shape A
$_{1}$ induces a larger phase lag, but in this case, the diminution of the tip amplitude is less important. Almost no phase lag is found for the shorter fins, independently of the trailing edge geometry (shapes A
$_{2}$ and C). The shorter geometries thus behave more like rigid fins. The capacity of each hydrofoil to generate streamwise fluid velocity can be analysed by averaging the value of 
$u_x$ over multiple periods and the whole volume (
$\overline {u}_{\mathit {vol}}$) or over the fins virtual boundary (
$\overline {u}_{\mathit {fin}}$), the latter option displaying more obvious differences (see table 2). The highest streamwise velocity is produced by the bilobed short fin (shape C), followed by the long narrow fin flapping at maximal rate (shape A
$_{1}$ at 
$f_{3}$). The slowest fluid velocity is observed in the case of the wider fin (shape B) followed by the long narrow fin flapping at minimal rate (shape A
$_{1}$ at 
$f_{1}$). Based on the Strouhal numbers alone, we would anticipate that the foil shape A
$_{1}$ with flapping frequencies 
$f_{2}$ and 
$f_{3}$ lie outside of the propulsive efficiency range [0.2, 0.4]. Nevertheless, the efficiency ratio 
$\eta$ offers a different perspective on that classification, as will be shown in § 3.2.
Hydrofoil midlines at five time points equidistant over half a period.
3.2. Spatio-temporal distributions of fluid forces
The colour code of figures 5, 6 and 9 follows that of figure 4 in terms of time partition. Figure 5 illustrates the distributions of fluid forces (per unit length) in the 
$x$ direction, acting on the fins, along their proximo-distal axis. Half a period is represented using five equidistant time points (naturally, the opposite half-period would produce symmetrical distributions). Negative values correspond to useful propulsive thrust as the system is pushed upstream. In all cases, the maximal amount of thrust is produced when the tip is leaving its maximal point and initiating its motion back towards the centre (time point 1/10 in green). The thrust then reduces gradually until time point 3/10 (orange), after which it increases again as the fin approaches the opposite extremity of its excursion (time point 4/10 in yellow). The maximal amount of thrust is generated by the long fin at maximal flapping rate (shape A
$_{1}$ at 
$f_{3}$), followed by that same geometry at frequency 
$f_{2}$. The wider fin (shape B) also produces significant thrust although slightly lower than the narrow geometry. For these three cases, the maximal amount of force is generated by the most distal portion of the fin, peaking close to 80 % of the total length. The thrust is drastically reduced in the case of the shorter fin (shape A
$_{2}$). In comparison, the short fin with bilobed edge (shape C) produces slightly larger amount of thrust at time points 0/10 and 1/10, but this is compensated by a positive distribution of 
$F_x$ at time 3/10, detrimental to propulsion. In that specific case, the maximum amount of force is generated closer to the centre of the foil. The time evolution of the total force in the 
$x$ direction (integrated over the whole fin surface) is shown in red in figure 7. Because of the symmetry between the left and right strokes of the fin, maximum thrust is generated twice per cycle.
Thrust (
$x$-component) force per unit length acting on the fin (in N m
$^{-1}$), at five time points equidistant over half a period.
Lateral (
$y$-component) force per unit length acting on the fin (in N m
$^{-1}$), at five time points equidistant over half a period.
Total forces (in mN) acting on the fin, at five time points equidistant over half a period (thrust and lateral forces, in the 
$x$ and 
$y$ directions respectively).
Figure 6 presents the 
$y$ component of the hydrodynamic force (per unit length) acting on the fin. This lateral force corresponds to wasted energy (not useful to propulsion). The most thrust-producing foils (shapes A
$_{1}$ at 
$f_{3}$, A
$_{1}$ at 
$f_{2}$ and shape B) are also the ones where the most energy in lost in accelerating fluid in the lateral direction. Moreover, the distribution on 
$F_y$ is correlated to the distribution of 
$F_x$ both temporally and spatially, with a maximum lateral force produced at time points 0/10, 1/10 and 2/10, at a position close to 80 % of the fin length. The time points which are typically not associated with thrust production (3/10 and 4/10) correspond to lower lateral forces. The bilobed geometry (shape C) presents a more complex spatio-temporal distribution of 
$F_y$, where the force peak alternates between two locations during the half period: at 50 % of the fin length for time points 1/10 and 3/10, and 80 % for time points 2/10 and 4/10. The lowest levels of wasted energy due to lateral forces are found for the long fin at low frequency (shape A
$_{1}$ at 
$f_{1}$), followed by the shorter fin (shape A
$_{2}$). The time evolution of the total lateral force is shown in blue in figure 7. The left and right strokes produce mirroring force profiles due to the motion symmetry.
As a measure of propulsive energy, we use 
$\eta$, the ratio between the energy employed to produce thrust and the sum of that energy with the energy employed to produce lateral forces (2.8). Figure 8 allows us to compare the efficiency thus defined for all foils with their respective flapping rates, as a function of 
$St_r$ and 
$St$. David et al. (Reference David, Govardhan and Arakeri2017) have raised the question of whether or not the Strouhal number calculated from the tip excursion amplitude is a good metric for the width of a vortex wake. They showed that for highly flexible foils, the width of the wake is overestimated by the large tip excursion, and that a rigid projection of the pitching peduncle approximates better the vortices spacing. Interestingly, the fact that flexible fins have larger tip excursion amplitudes than their rigid counterparts is not verified in all our experiments. It is true for geometries A
$_{1}$ and A
$_{2}$, but not for shapes B and C, where the amplitude is smaller than the rigid projection. Therefore, in some cases, the flexibility of a fin tends to decrease its excursion amplitude and its Strouhal number.
Efficiency ratio (as defined by (2.8)), as a function of the rigid Strouhal number 
$St_r$, calculated based on the rigid projection of the peduncle (a) and the flexible Strouhal number 
$St$, calculated based on the tip excursion amplitude of the flexible foils (b).
The most efficient shape is the wider foil (B), followed by the long narrow shape (A
$_{1}$), which is equally efficient at frequencies 
$f_{2}$ and 
$f_{3}$ but drops in efficiency at lower frequency (
$\,f_{1}$). The shorter fins are the less efficient ones, with the lowest value of 
$\eta$ found for the short fin with a bilobed edge (C). This is consistent with the fact that flexible fins are typically more efficient than their rigid counterparts (Dewey et al. Reference Dewey, Boschitsch, Moored, Stone and Smits2013). Indeed, shapes A
$_{2}$ and C, because of their specific geometries, behave almost like rigid panels, as can be inferred from the close to zero phase lag between the peduncle and the tip angles (table2), and from the visualization of their midlines (figure 4).
In figure 8, the dashed regions denote the range of Strouhal numbers typically associated with efficient thrust production. Interestingly, the classification based on the definition of 
$\eta$ shows that shape A
$_{1}$ at frequencies 
$f_{2}$ and 
$f_{3}$ is more efficient than shapes A
$_{2}$ and C, even though their Strouhal numbers lie outside of the usual window of efficiency [0.2,0.4]. In all cases, the efficiency ratio 
$\eta$ remains below 0.2, indicating that less than 20 % of the input power is used for propelling the system in the upstream direction, which is not particularly efficient.
The normal force distributions (per unit length) on both sides the hydrofoils are illustrated in figure 9. A low flapping rate (shape A
$_{1}$ at 
$f_{1}$) or a shorter fin length (shape A
$_{2}$) results in lower normal forces across both surfaces, peaking close to 70 % of the total fin length. At intermediate or higher flapping rate (shape A
$_{1}$ at 
$f_{2}$ or 
$f_{3}$), when the fin is wider (shape B), or when the edge is bilobed (shape C), the normal forces increase and the maximum magnitude is shifted close to 80 % of the fin. Moreover, the asymmetry between the left and right sides of the hydrofoils becomes more important. The magnitudes are asymmetric, and in the most extreme cases, the normal force has the same sign on both fin sides (see for example time points 3/10 and 4/10 for shape A
$_{1}$ at 
$f_{2}$). These asymmetries translate into larger internal tension within the material. If the sum of the forces on both sides is positive, the material inside the fin is subjected to an outward stretching as it experiences a net force in the direction of the outward normal. Contrarily, if the sum on both sides is negative, the resultant force points towards the inside of the fin, which is then subjected to a net compression. A perfectly symmetric force distribution on both sides (sum equal to zero) would indicate that one side is pulled while the other is pushed with the exact same magnitude, therefore, the hydrodynamic stress would result in the acceleration of the foil, with no internal stress caused by the net external contact forces. To evaluate the spatial distribution of internal tension across the hydrofoils, the normal stress is averaged over a full cycle of flapping, with its sign. The results are presented as colour maps in figure 10. Negative values denote an excess of compression over the oscillation cycle, positive values denote stretching. All membranes experience an overall stretching over the tip region. The most extreme occurrence of that effect is found at the tip of shape A
$_{1}$ at 
$f_{3}$. A high concentration of stretching is also seen for shape A
$_{1}$ at 
$f_{2}$, localized in the very centre of the tip, and for shape B, where the concentration zone is shifted to the corners. Additionally, there is a compression zone located distally to the centre of the foils, between approximately 50 % and 75 % of the total length. This is seen in all cases except for shape A
$_{1}$ at 
$f_{2}$, where the internal compression is relaxed in the centre of the fin. Furthermore, the normal stress distributions are averaged over a period, in absolute values, revealing which portions of the fin are most prominently involved in the production of hydrodynamic forces. These spatial distributions are presented in figure 11, with stress curves covering the dorso-ventral axes of the fin and taken at different locations across its length (25 %, 50 %, 75 % and 100 %). The maximum values are always found at 
$\sim$75 % of the foil length. The variations across the surface are marked mostly in the case of the long and narrow geometry (shape A
$_{1}$), where the curves at 0.25
$L$, 0.5
$L$, 0.75
$L$ and 
$L$ become more and more separated as the flapping frequency increases. It is noteworthy that a special stress pattern emerges for shape A
$_{1}$ at 
$f_{2}$, where the curve at 0.75
$L$ adopts a bilobed signature, not observed for 
$f_{1}$ nor 
$f_{3}$. The excess of normal stress on the tip corners of shape B is also visible in the corresponding graph (stress curve at 
$L$). Finally, all fin geometries and flapping rates result in lower levels of hydrodynamic forces generated at the lateral edges (dorsal and ventral) of the hydrofoils.
Normal force per unit length (in N m
$^{-1}$), acting on both sides of the fin (left side: diamonds, right side: stars) at five time points equidistant over half a period.
Period-averaged signed value of the normal stress acting on the sides of the fin (in Pa).
Period-averaged absolute value of the normal stress acting on the sides of the fin (in Pa, 
$\sigma _p=0.2$ Pa).
3.3. Vorticity fields and wake structures
The vorticity field (
$\boldsymbol {\omega }=\boldsymbol {\nabla } \times \boldsymbol {u}$) is very informative about characteristic flow structures, revealing relations between wake topology, surface distributions of fluid forces and propulsive efficiency. To allow for comparison with previous work, we present the distributions of 
$\omega _z$, the 
$z$ component of the vorticity, in 2-D planes intersecting the fins at mid-span. Instantaneous 
$\omega _z$ maps are shown in figure 12, as the foils are crossing the 
$y=0$ axis (corresponding to 2/10 of a period, colour coded in pale blue in figures 4–6 and9). Period-averaged 
$\omega _z$ maps are shown in figure 13.
Snapshots of 
$\omega _z$, the 
$z$-component of the vorticity field (in s
$^{-1}$) in the plane intersecting the fins half-way through their span, as they are crossing the 
$y=0$ axis, travelling towards positive 
$y$.
Period average of 
$\omega _z$, the 
$z$-component of the vorticity field (in s
$^{-1}$) in the plane intersecting the fins half-way through their span, with the cycle of motion envelopes indicated by white areas.
During the flapping cycle, the instant when the fins cross the 
$y=0$ axis occurs between the points of maximal and minimal thrust production (see figure 7). Figure 12 illustrates that a vortex with a core of negative 
$\omega _z$ (rotating clockwise) forms at the tip of each fin. For shape A
$_{1}$ at frequency 
$f_{1}$, the vortex is almost detached at this point, whereas at 
$f_{2}$, it is fully detached, with a second smaller vortex forming at the tip. The size and strength of these clockwise vortices correlate with propulsive efficiency: the largest vortices, with maximal 
$| {\omega _z}|$ of the order of 50 s
$^{-1}$, are found for shape A
$_{1}$ at 
$f_{2}$, shape A
$_{1}$ at 
$f_{3}$ and shape B at 
$f_{2}$, the three cases with maximal 
$\eta$ coefficients. Counterclockwise vortices are also present in the snapshot vorticity fields for shapes A
$_{1}$ (at 
$f_{3}$), A
$_{2}$, B and C. They are located more downstream as they have been shed during the precedent half-cycle of motion. In the lower frequency cases for shape A
$_{1}$ (
$f_{1}$ and 
$f_{2}$), they are not discernible in the field of view, most likely because they have been shed downstream past 
$x=35$ mm by that time. The persistence in the close downstream wake of these positive 
$\omega _z$ vortices do not appear to be correlated with propulsive efficiency, nor their magnitude or orientation. Rather, it is the flow structures directly adjacent to the tip which embody the capacity of a fin to generate efficient thrust. It is interesting to observe that the curvature of the trailing edge has an impact on the position at which the vortices are being shed. Indeed, counter-rotating vortices are added to the wake along the 
$y=0$ axis for a straight edge geometry (shape A
$_{2}$), whereas they are detaching at higher 
$y$ positions for a bilobed trailing edge (shape C).
The period-averaged 
$z$-vorticity also presents a correlation with propulsive efficiency: the three hydrofoils with highest 
$\eta$ coefficients (shape A
$_{1}$ at 
$f_{2}$ and 
$f_{3}$ and shape B at 
$f_{2}$) display the highest means, approximately 15 s
$^{-1}$ (figure 13). For all shapes and frequencies, the period-averaged 
$\omega _z$ indicates that a surplus of clockwise vorticity is induced at the tip when the fin travels in the 
$y<0$ region, and vice versa. In other words, the fluid tends to rotate more towards the exterior of the motion envelope close to the tip. In the wake surrounding more proximal portions of the foils, the opposite effect is observed: the period-averaged 
$z$-vorticity is slightly positive in the 
$y<0$ region, and vice versa. This phenomenon is linked to the period-averaged surface maps of normal stress (figure 10), denoting an excess of stretching close to the tip, and a surplus of compression towards the centre and more proximal surface regions. Moreover, a stronger stretching force at the tip appears to be related to higher mean vorticity close to the trailing edge, as seen for shape A
$_{1}$ at 
$f_{2}$ and 
$f_{3}$ and shape B at 
$f_{2}$. From figure 13, an estimation of the vortex shedding angle can be drawn, corresponding to the opening angle of the maximal 
$ |{\omega _z}|$ zone. This angle is wider in the cases of shape A
$_{1}$ at 
$f_{2}$ and 
$f_{3}$, roughly following the tip curvature of the motion envelopes. The most efficient fin (shape B) presents a less pronounced shedding angle, which remains aligned with the motion envelope outline. The shedding angle seems particularly constrained for the two shorter fins behaving more like rigid panels (shapes A
$_{2}$ and C), which are also the least thrust-efficient hydrofoils. In those cases, the high vorticity areas are almost aligned with the 
$y=0$ axis. When comparing the period-averaged wake signatures of shape A
$_{1}$ at 
$f_{1}$, 
$f_{2}$ and 
$f_{3}$, we observe that increasing the flapping frequency results in the appearance of a reversed circulation area in the downstream wake, close to 
$x=33$ mm, where 
$\omega _z$ changes its sign. This mean vorticity reversal does not seem related to propulsive efficiency, since it is also slightly visible in the case of shape A
$_{2}$ at 
$f_{2}$, one of the least efficient foils in terms of thrust production.
4. Discussion and outlook
We have established a systematic experimental approach to compare the hydrodynamic stress maps on different fin shapes with varying flapping frequency, allowing us to identify the effects of each parameter on the spatio-temporal distributions of fluid forces, and link them to the direct wake topology. In our analysis, we payed special attention to the ability of each hydrofoil to generate thrust (force propelling the system upstream), lateral forces (not contributing to propulsion) and the resulting propulsion efficiency. We also investigated the spatial distributions of external contact forces normal to the fins surfaces and the corresponding internal tension accumulating within the material. We have found that the long narrow geometry generates the most thrust (at equal frequency), and that increasing the flapping frequency increases thrust. Between both short geometries, the bilobed trailing edge allows us to produce slightly more thrust than the straight edge, but only temporarily during the period. At a certain point during the flapping cycle, it is also pushed downstream by the fluid (the opposite of propulsive thrust). It is typically the distal region of the fins which participates more in thrust generation, except for the short bilobed geometry whose central portion is the one most involved. As for the lateral forces, detrimental to propulsive efficiency, the short straight geometry (at equal frequency) reaches the lowest levels. For the long geometry, reducing the flapping frequency also reduces the lateral force magnitude throughout the cycle.
The examination of instantaneous and period-averaged vorticity fields yielded the conclusion that the formation of larger and stronger vortices directly at the tip as well as the accumulation of excess vorticity close to the trailing edge over a full period are related to better propulsive efficiency. Among all cases tested, the most efficient geometry was found to be the long wide fin, followed by the long narrow one, whose efficiency remained similar at higher frequency but dropped significantly at the lowest one. The least efficient fin morphology is the short fin with bilobed trailing edge. This type of tip curvature is ubiquitous in nature. It is interesting, both for the study of biological swimmers and for the design of artificial fins, to discover that the naturally widespread bilobed shape does not improve efficiency nor does it grant more thrust producing capabilities to a simple fin model at Reynolds numbers which are on the low side of the inertial regime. It is more likely that an authentic fish fin with a bilobed trailing edge is coupled with specific elasticity profiles and complex controlled motion programs, in order to achieve its functional role. It could also be that the advantages of this tip curvature appear only for a specific combination of parameters such as length and flapping frequency. Besides, we have observed that the bilobed geometry results in a dynamical shift of fluid force peaks along the fin length during the stroke cycle. This more complex dynamics of force distributions could be used actively by an appendage to achieve fine manoeuvres. The external contact forces from the fluid acting normally on both foils surfaces result in an accumulation of internal tension throughout the flapping cycle. The fin tip is subjected to an excess of tissue stretching, an effect which becomes more important as the foil becomes longer, or as it flaps faster. The stretching tension concentrates towards the tip corners for a wider geometry. Moreover, the foil membranes experience an overall compression zone at approximately two thirds of their total length. Interestingly, in the case of the long narrow geometry, flapping at frequency 
$f_{2}$, this compression region partly disappears. In that specific case, a particular stress pattern emerges from 75 % of the total length up to the tip, where the normal stress (averaged in absolute value over a period) presents a double peak signature. For all fin geometries, the period-averaged absolute value of the normal stress reaches its maximum at approximately 75 % of the total length. These stress maps are very indicative of how the fin architecture must be conceived in order to withstand specific concentrations of mechanical loads and internal tension throughout its oscillation cycle.
The flexible tails of rayed-fin fishes display an impressively wide range of morphological traits and provide a great model system to study the relation between fin geometry and propulsive efficiency (Affleck Reference Affleck1950; McNeill Reference McNeill1974; Lauder Reference Lauder1989; Weihs Reference Weihs1989; Lauder Reference Lauder2000; Lauder & Madden Reference Lauder and Madden2007; Blake, Li & Chan Reference Blake, Li and Chan2009; Flammang et al. Reference Flammang, Lauder, Troolin and Strand2011b). Complex combinations of environmental constraints participate in the selective evolution of fish fins, and the resulting elastic properties, geometry and kinematics of each specie constitute an elegant compromise between competing survival demands such as the necessity to achieve fine control of the fin surface for delicate manoeuvring, to generate quick acceleration bursts or to maintain long distance efficiency (Webb Reference Webb1975; Lauder & Liem Reference Lauder and Liem1983; Lauder Reference Lauder1989, Reference Lauder2006, Reference Lauder2015). It is tempting to assume that evolution has optimized the geometric and kinematic parameters of fins to meet with the functional demands of distinct habitats. However, drawing a parallel between form and function requires taking into account the evolutionary history of fins morphologies, and considering not only the external shape, but also internal structures, in order to avoid mistaking correlations between morphology, habitat and function for causal relationships (Lauder & Liem Reference Lauder and Liem1983; Lauder Reference Lauder1989). Previous investigations on the interplay between shape and stiffness of fin models and their capacity to generate thrust efficiently have already demonstrated that efficiency varies in a complex manner as a function of the different fin attributes, and that it is not a simple matter to identify an optimal shape (Feilich & Lauder Reference Feilich and Lauder2015). Rather, efficiency can be maximized by several combinations of parameters, depending on the flow conditions, and it remains a great challenge to draw general conclusions about the functional advantages of geometrical and kinematic aspects. The present work illustrates how the investigation of hydrodynamic stresses based on three-dimensional PTV experiments can provide precious information about the mechanical interplay between a fin-like structure and the surrounding fluid. It paves the way for future experimental studies using direct evaluation of fluid forces on the surface of submerged flexible appendages, allowing us to tackle directly the fundamental problem of form and function in the context of propulsion through a fluid medium. To conclude, the possibility of resolving the 3-D instantaneous hydrodynamic stress fields on small-scale models of deformable structures opens the door to many applications in the engineering design of aeronautic and aquatic vehicles or the field of biophysics concerned with animal locomotion.
Acknowledgements
This work was funded by the Swiss National Science Foundation (SNF) via a Sinergia research grant as well as a UZH Forschungskredit Candoc grant. We are grateful for interdisciplinary discussions with T. Aegerter, A. Jaźwińska, I. Kicic, S. Puri and S. Verma. We are very thankful to J. O. Dabiri et al. for making the queen2 pressure algorithm available (Dabiri et al. Reference Dabiri, Bose, Gemmell, Colin and Costello2014).
Declaration of interests
The authors report no conflict of interest.
