Hydrodynamic stress maps on the surface of a flexible fin-like foil

We determine the time dependence of pressure and shear stress distributions on the surface of a pitching and deforming hydrofoil from measurements of the three dimensional flow field. Period-averaged stress maps are obtained both in the presence and absence of steady flow around the foil. The velocity vector field is determined via volumetric three-component particle tracking velocimetry and subsequently inserted into the Navier-Stokes equation to calculate the total hydrodynamic stress tensor. In addition, we also present a careful error analysis of such measurements, showing that local evaluations of stress distributions are possible. The consistency of the force time-dependence is verified using a control volume analysis. The flapping foil used in the experiments is designed to allow comparison with a small trapezoidal fish fin, in terms of the scaling laws that govern the oscillatory flow regime. As a complementary approach, unsteady Euler-Bernoulli beam theory is employed to derive instantaneous transversal force distributions on the flexible hydrofoil from its deflection and the results are compared to the spatial distributions of hydrodynamic stresses obtained from the fluid velocity field.


Introduction
The evaluation of hydrodynamic forces on the surface of submerged and deforming foils constitutes a key element in understanding how synthetic and animal fins interact with the surrounding fluid to achieve their appropriate functions of propulsion and maneuvering. Great effort is deployed by scientists to describe the motion of fish fins and explore the physics underlying their complex kinematics [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. An experimental workflow is proposed here to measure local distributions of hydrodynamic pressure and viscous stress over the surfaces of a flexible artificial fin, starting from the volumetric acquisition of 3D flow velocimetry data. This methodology broadens the possibility of studying in detail the propulsion and force regulating mechanisms in fish or biomimetic systems [18].
Besides, several groups have resorted to PIV/PTV methods to resolve the laminar and turbulent shear layers forming over solid surfaces moving in water, such as undulating membranes [94] and fish [108]. Wall shear stress determination based on PIV was also conducted in small arteries models (diameter below 1 cm) with blood mimicking fluid [109]. Two important difficulties arise in the evaluation of wall shear stress from the fluid velocity field: (1) the ratio of boundary layer thickness to spatial resolution must be large enough to allow the estimation of the velocity gradients inside that layer and (2) the stress profiles must be measured as close as possible to the real surface, where reflections can hinder the detection of particles. The shear stress profile remains fairly well captured (with a slight offset in magnitude) as long as the reconstructed surface stands sufficiently close to the real object (distance below 1 mm according to [110]). Using volumetric PTV experiments, we recently showed the feasibility of resolving the velocity field over a rigid plate pitching inside a water tunnel, with a spatial resolution of 0.75 mm, and we measured the thickness of the boundary layer within a dynamic range of about [2,6] mm [111].
In the present work, we apply a similar methodology to resolve the shear forces over the surface of a flexible hydrofoil flapping in comparable flow regimes, combined with the pressure gradient integration method. Thus, we demonstrate the possibility to obtain well resolved distributions of pressure and viscous shear stresses, both spatially and temporally, on the surface of a synthetic fin whose size is of the order of 1 cm 2 . For this purpose, we apply a full volumetric, three components particle tracking velocimetry technique (3D-3C PTV) for the first time in this context. In our study, special attention is paid to the time and spatial resolutions and their experimental limitations, and we present a careful analysis of the error propagation from the tracked particle positions to the final hydrodynamic stress values. This work demonstrates the feasibility of experimentally capturing the full 3D force field in the fluid environment of a flapping flexible foil with dimensions comparable to that of a small fish, and to employ a force decomposition method to project the hydrodynamic stress tensor on the relevant morphological axes of the fin.
As a complementary calculation, we use the unsteady Euler-Bernoulli beam deflection theory to derive the instantaneous transversal force distributions acting on the deflecting hydrofoil. This approach relies on the assumptions of small angle deflections and rigid, planar crosssections, remaining perpendicular to the deformed midline, considering deflections in one plane only [112]. Therefore, it offers a more approximative framework to describe fluid forces acting along the foil, as compared to the hydrodynamic forces derived from the 3D flow field. Nevertheless, both methods have been employed in previous work to study flexible foils flapping inside a fluid, and it is interesting to compare their results for the same experiment.
The rest of the paper is organized as follows. In section 2.1, the basic principles behind the experimental setup are described, before the characteristics of the model fin and the hydrodynamic stress calculation are introduced in the subsequent sections. In order to calculate forces on the respective fin surfaces from the stress tensor, the surface needs to be parameterized, which is described in section 2.4. The description is continued with the basis of alternative methods, such as Euler-Bernoulli theory and a control volume analysis, before the materials and methods section is closed with a description of the limitations in uncertainty and resolutions. Section 3.1 shows the main results in the determination of the hydrodynamic stress tensor for two different experimental situations, which are validated in a control volume analysis in section 3.2 and compared to results from Euler-Bernoulli theory in section 3.3. Finally, the limits of the method are presented in terms of the uncertainties for the obtained forces and how they relate to different flow structures in section 3.4. The results are put into context with possible applications in the discussion and outlook.

Volumetric particle tracking velocimetry
PIV and PTV techniques have become widespread in the past decades in various hydrodynamics research fields and have been extensively described in the literature [113][114][115][116]. We performed volumetric three-components PTV measurements using the V3V-9800 system (TSI Incorporated, 500 Cardigan Road, Shoreview, Minnesota 55126 USA), which is well characterized in [117]. This method presents the advantage of capturing the full unsteady flow field in a single recording of the whole 3D domain, and thus is preferable to more arduous implementations of scanning stereo-PIV [63, 95,118], where a large number of measurement planes is required to reconstruct the whole volumetric flow field, or tomographic PIV, where the whole volume needs to be scanned with a light sheet for each particle position field, which is challenging with a flexible membrane deforming rapidly in the flow [66,119].
In the present setup, illustrated in Fig 1, a synthetic fin is attached to a small metallic rod and inserted through the back wall of the flow chamber, with a servomotor fixed outside the water tunnel to activate the pitching program of the rotation axis. Double pulses emitted by a dual laser system illuminate a measurement volume (50×50×20 mm 3 ) inside the flow chamber. The working fluid is water and the seeding tracers are polyimide particles with a diameter of *50 μm, which were used in a precedent study of fish locomotion based on particle velocimetry [15]. In total, *45 mL of particles are seeded in *160 L of water, which yields an approximate number of 8×10 4 particles inside the interrogation volume, considering that the spherical particles powder has a filling fraction of *40%. Three cameras mounted on a triangular plate are recording triplets of images. Hence, the tracer particles advected by the flow are captured from three different perspectives at each time point. Each of the three images is analyzed based on a 2D Gaussian fit of the particle intensity distribution to identify the 2D position map. The 3D field of particle positions is reconstructed by overlapping the three viewing angles based on a triplet search algorithm. Flow velocities are determined based on individual particle displacements in the time between subsequent laser pulses (δt = 2.5 ms), using a relaxation method to achieve a probability-based matching. Appropriate synchronization is arranged between the cameras and lasers in a temporal scheme called frame straddling, as illustrated in Fig 1. The rate at which velocity fields are recorded (80 Hz) is half the camera acquisition frequency, yielding a time separation of Δt = 0.0125 s between the velocity fields. This allows the reconstruction of the material acceleration field, necessary for the hydrodynamic stress calculation.
The image post-processing was conducted using the V3V software (version 2.0.2.7) along with an optimal set of parameters, selected to reduce the number of velocity outliers while preserving the main flow features. The maximum overlap between two particle intensity distributions was set to 65%, a mask was applied over the foil during the particle identification step, and appropriate median filter and velocity range (±0.25 m/s) were applied to the raw velocity (1) Flow chamber with transparent windows on three sides and a fixation wall on one side for inserting the synthetic fin. (2) Dual-head pulsed Nd:YAG laser with wave-length of 532 nm and maximal energy of 120 mJ/pulse, equipped with a pair of cylindrical lenses to expand the beam. (3) Mirror to deflect the laser beam. (4) Three cameras (resolution 4 MP, 85 mm lenses, sensor size 11.3×11.3 mm 2 , magnification 0.3, maximal frequency of capture 180 Hz) mounted on a triangular plate parallel to the front wall of the flow chamber (x-y plane). The distance between the cameras plate and the center of the water tunnel is *465 mm. (5) Water tanks for the recirculating system. (6) Pipe and pump to control the flow. (7) V3V software and synchronizer. (8) Close-up of the flow chamber with inner dimensions (in cm), PTV-interrogation volume in blue (*50×50×20 mm 3 ) and sketch of the trapezoidal fin. Frontal view (x−y plane) and top view (x−z plane). (9) Frame straddling method with δt between the particle images and Δt between the velocity fields. https://doi.org/10.1371/journal.pone.0244674.g001 fields. Finally, the particle velocities were interpolated on a regular 3D grid using Gaussian interpolation, yielding a final spatial grid spacing of 0.75 mm in each direction. Precedent studies found that appropriate smoothing of the velocity fields could reduce the error in pressure fields produced by the queen2 algorithm from [77] by 30-67% [80,107]. Precautions need to be taken for turbulent flow fields to avoid smoothing out high frequency flow structures which could compromise the pressure field reconstruction. In this case, a moderate Gaussian smoothing scheme with a factor of 1.5 (corresponding to the ratio between the Gaussian radius and the voxel half size) was applied to the interpolated velocity field. 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 (see section 2.4).  Table 1. The flexible foil was produced by 3D printing a rigid cast of its negative form, then pouring liquid PDMS (polydimethylsiloxane) inside the cast, installing the fork-shaped base rod at the appropriate location, and curing the model for 36 hours at 58˚C. This resulted in a flexible membrane of thickness 0.55 mm. The mass density fo the foil ρ f matches that of water (http://www.mit.edu/~6.777/matprops/pdms.htm). A study of elastic properties of real caudal fins has reported an effective Young's modulus of *8 MPa for the zebrafish, including both the rays and interray tissue contributions [120,121]. The same cantilever deflection setup was employed by Sahil Puri (University of Zurich) to characterize the  elastic properties of our PDMS foil, resulting in a value of *0.8 MPa. In order to further characterize the properties of our fin model, we performed a mechanical damping analysis, where the fin was simply held from the tip at a certain amplitude and released (in air and in water). The time evolution of the tip amplitude was modeled by a decaying oscillation:

Geometry and kinematics of the synthetic fin
We found values of γ air = 1.5 s -1 and γ water = 5.0 s -1 for the decay rates in air and in water, respectively, and values of f 0,air = 9.3 Hz and f 0,water = 2.4 Hz for the eigenfrequencies. These results are needed for the estimation of the damping terms in the Euler-Bernoulli equation (see section 2.5).
Two versions of the pitching fin experiment were conducted, with different pitching amplitude θ 0 and upstream velocity u 1 ( Table 1). The dimensionless Reynolds and Strouhal numbers are employed to characterize the flow regime of the flapping foil: These parameters are based on the kinematic viscosity of water (ν), the fluid velocity averaged over the whole measurement volume and a complete flapping cycle (ũ ave ), the foil length (L = 25 mm), the flapping frequency (f = 3 Hz) and the tip-to-tip excursion amplitude (a, obtained from the parametrized deflection profiles-see section 2.4). Even within a single fish species, a wide range of St and Re can be reached by caudal fin motion, depending on the developmental stage (larval, juvenile or adult) and the swimming behavior. For example, zebrafish undergoes a developmental transition which allows it to cover a large spectrum of Reynolds and Strouhal numbers throughout its lifetime, ranging from 390 to 2600 in Re (using the caudal fin length to allow direct comparison with isolated hydrofoils) and from 0.3 to 2.5 in St, based on swimming speed, tail beat frequency and amplitude reported by numerous studies [14,17,[122][123][124][125][126][127][128][129]. The Reynolds numbers of the two experimental cases presented here (385 and 725), although not representative of maximal velocity in fish locomotion, correspond to slow cruising regimes characteristic of certain species [130] and specific behaviors such as foraging, chemotaxis and exploration. The chosen flapping frequency (3 Hz) is also comparable to typical tail-beat frequency in several fish species [131][132][133][134]. Moreover, the frequency was selected in the vicinity of the resonance frequency of the foil-fluid system (f 0 = 2.4 Hz). In addition to the well known Strouhal number efficiency window of (0.2-0.4) [135,136], it is generally accepted that a driving frequency close to the eigenfrequency of the oscillator will result in improved efficiency in thrust production [38,65,[137][138][139][140][141]. Therefore, experiments with flexible foils flapping close to the resonance frequency are of particular relevance for the field of underwater propulsion and fish locomotion. The frequency chosen for our experiments correspond to a flexural frequency of 1.25 (defined as f � = f/f 0 ). The flexural rigidity f � characterizes the relative time scale of the flapping motion with respect to the time scale of the fluid-structure resonance. It can be used to describe the fluid-structure interaction, similarly to the Strouhal number which characterizes the relative time scales of flapping behavior with respect to the external flow [141]. Therefore, although not yielding values of St typically related to high propulsion efficiency, the chosen value of the flapping frequency is pertinent for hydrodynamic problems concerned with the behavior of flexible structures oscillating near resonance, and is comparable to values employed in other studies of fin-like hydrofoils [80,137,139,[141][142][143][144][145][146]. The experimental cases presented here are not designed to achieve high propulsive efficiency, rather, to demonstrate the feasibility of obtaining well resolved stress maps on all fin surfaces. This methodology could be applied to other ranges of St and Re to study problems concerned with specific flow regimes, for example, by tuning the external incoming flow velocity or the flapping frequency, as shown in a recent study [18].

Hydrodynamic stress calculation
The starting point of the hydrodynamic stress calculation is the Navier-Stokes or momentum equation [147,148]: The last term on the right corresponds to any type of body force acting on the fluid, such as gravity which can be measured independently of the flow. This can be included in a static pressure term and we will omit it in the following. The left side of the equation contains the material acceleration, with the Lagrangian derivative D/Dt. The first term on the right side is the divergence of the total hydrodynamic stress tensor s ij = −pδ ij + τ ij , where p is the scalar pressure field and τ the viscous stress tensor. This viscous shear stress can be expressed as: where μ is the dynamic viscosity. Computationally, τ ij can be obtained directly using centered finite differences between neighboring velocity vectors: In order to determine the pressure from the flow fields, we first rewrite the Navier-Stokes equation to isolate the pressure gradient on one side: The pressure field can then be reconstructed by integrating this equation, with the additional requirement of specifying the pressure values at the starting points of the integration paths pðr ref Þ.
The pressure calculation was conducted using the queen2 algorithm [77], which is available at http://dabirilab.com/software/. The calculation assumes zero-pressure values on the external boundary of the domain, in the undisturbed flow, and integrates the pressure gradient along eight different paths (horizontal, vertical or diagonal) originating from the outside contour, among which the median is selected for each node in the 3D domain. The object is masked to prevent any integration path crossing its boundary. This algorithm offers the advantage of reasonable computation times even for large domains, with a total number of paths equal to 8 times the domain size (N x × N y × N z ) per velocity field, resulting in a more than order-of-magnitude improvement in computation time compared to other multi-path integration schemes.
Omni-directional schemes can provide higher accuracy [69,75,84,105] but for 3D implementation, the use of GPU-based calculation was demonstrated to achieve accelerated computational times [106]. The terms containing out-of-plane derivatives of the velocity (grouped separately in Eq 6) are not included directly in the computation. Rather, the three-dimensionality is handled through a combination of integration paths originating from all the faces of the domain, combined with a median polling. More details on how this median polling approach allows to reduce the accumulation of measurement errors are provided in [77].
The central step of the pressure calculation consists in evaluating the material acceleration, which can be done in the Eulerian or the Lagrangian frame. The Eulerian frame is stationary in the laboratory, where the time derivative of the velocity is expressed as a total derivative: In the queen2 algorithm, the material acceleration is evaluated based on a Lagrangian forward scheme (Fig 3), which consists in following a particle trajectory in a pseudo-tracking formulation: The particle velocityũ � p ðx � p ðt iþ1 Þ; t iþ1 Þ at a forward instant relies on a linear reconstruction of the particle path to determine its positionx � p ðt iþ1 Þ, based on a trapezoidal scheme: An improvement of the pseudo-tracking approach was proposed in [69], relying on two cameras to obtain combinations of four staggered exposures of the particle traces in order to determine the Lagrangian velocity. The simplified Lagrangian scheme from [77] is used instead, which was shown to provide accurate time-resolved pressure contours for different numerical and experimental cases, including a simulated 3D anguilliform swimmer [77], and most relevant for the present study, a flexible foil flapping inside a water tunnel, where the total forces on the foil compared favorably with load cells measurements [80].
Two types of stress maps are extracted from the flow field: instantaneous distributions and period-averaged distributions, i.e. averaged over equidistant time points across a full motion cycle. In experiment 1, eight pairs of velocity fields covering a complete period are selected for the calculation of instantaneous stresses. As a second step, the stress tensor is projected on the fin surfaces, and the forces are decomposed into biologically relevant components (along the fin principal axes-see Fig 2). The normal stress (in absolute value) is averaged for these eight time points to obtain the period-averaged hydrodynamic load over the fin. In experiment 2, the period-averaged normal stress (in absolute value) is presented to offer a comparison with the results from experiment 1. This allows to observe which portions of the fin surface are subjected to the highest normal loads, in average during a full cycle, and to determine how this distribution changes with or without an external current (experiment 1 versus experiment 2).

Surface parametrization
The raw images offer a very clear view of the foil midline curvature (see left panel of Fig 4 for example), and are used to parametrize the midline position with a polynomial of degree 2. For each instantaneous time frame under study, 25 points are manually superimposed on the fin, allowing to fit the quadratic function to the hydrofoil time-dependent deflection. The 3D hydrofoil, including the left and right trapezoidal sides plus the tip surface, is reconstructed based on the parametrized midline, assuming that deflection occurs only along one axis. This is a reasonable assumption since the foil length is about twice as large as its width. The virtual surface reconstructed from the parameterized deflection is located at a small distance from the real foil (0.725 mm) owing to the fact that the particles can not be resolved directly at the surface due to reflections (see the masked region around the fin in Figs 4 and 7 for example).
A crucial step in our analysis is the projection of the total hydrodynamic stress tensor onto the different surfaces of the foil. The i-component of the stress (force per unit area) acting on a surface with a unit normal vectorn (oriented outwards) is expressed as: The stress vector is initially expressed in cartesian coordinates. It is then decomposed into more relevant components, related to a set of orthogonal axes moving along with each fin surface, as illustrated in Fig 2. Anticipating useful analogies between our fin model and a real caudal fin, we employ the following terminology to define this coordinate system: the dorsoventral axisn dv , the proximo-distal axisn pd and the left-right axisn lr . These orthogonal sets are defined locally everywhere on the foil surfaces (sides and tip) and travel with it during the motion cycle. Thus, each grid point on the foil surfaces (25 nodes along the dorso-ventral axis, 25 nodes along the proximo-distal axis, 3 nodes along the left-right axis) is associated with a local set of perpendicular unit vectors (n pd ,n lr ,n dv ), on which the stress vector is projected to give S pd , S lr and S dv respectively. The pressure contribution to the total stress is always automatically acting along the normal direction onto the local surface.
In the present experimental cases, the fin operates in the inertial flow regime (although theoretically very close to the transitional range 300<Re<1000 [149]), where the normal stress component (dominated by the pressure) is typically larger than the viscous tangential stresses by at least two orders of magnitude. Depending on the flow under investigation, these distinct stress components may have different relative magnitudes. For example, the fish larvae evolve in a flow regime where the viscous shear stress is substantially larger [149]. Moreover, the tangential shear stress may have a particular significance in certain hydrodynamic problems, even when the flow is dominated by the fluid inertial forces. For instance, resolving the viscous shear stress maps on the surface of fin models might be informative in problems related to flow sensing in the skin of fish [111]. In other words, despite the lower order of magnitude of the measured viscous tangential stress, these components may be of special interest for certain biological applications, which is why we retained them in our analysis.
Vast effort has been dedicated in refining PIV/PTV techniques to research complex flow phenomena in boundary layers, including recirculating zones in cavities [69], cavitation [71], pressure fluctuations [78,98] and attached/separated layers over obstacles [150]. In the present work, the focus is on providing an estimate of both hydrodynamic stress contributions (shear stress and pressure) from instantaneous captures of the three-dimensional flow field, which requires a compromise between the size of the measurement volume and the vector field spatial resolution. In a recent study, we showed that it is possible to resolve the velocity field inside the shear layer forming over a flat plate pitching in a similar flow regime using the same volumetric PTV experimental setup [111]. In that case, the boundary layer thickness, defined as the distance at which the velocity reaches 99% of the free-stream value, varied dynamically between *2 and 6 mm. In the case of our flexible hydrofoil, a rough estimate of the boundary layer thickness can be calculated from the solution for a laminar layer over a flat plate [151]: The value of x is the distance along the proximo-distal axis measured from the foil leading edge. The stress maps and curves presented in section 3.1 correspond to locations x on the fin in the range of [0.25L, L], where L is the total length of the foil. Even for the thinnest boundary layer, which should occur in exp. 2 where the average velocity magnitude is the largest, the boundary layer thickness evaluated by Eq 11 lies in the range of [2.3,4.6] mm, which is at least three times larger than the vector field spatial resolution (0.75 mm). Despite potential transitions to turbulence due to the oscillation of the flexible membrane, it is reasonable to assume that this thickness will remain in a similar range as that measured on the pitching flat plate [111], thus allowing the estimation of velocity gradients inside the shear layer. Additionally, the virtual surface distance (0.725 mm) is small enough compared to the boundary layer thickness so that the measured tangential stresses can reflect the trends inside the shear layer [110].
This empirical approach allows the estimation of all hydrodynamic force contributions from a single time series of three-dimensional flow fields, thus offering a relevant approach for problems requiring simultaneous shear stress and pressure determination, for example, on the surface of a fish body or fin whose boundary layer thickness is of similar order of magnitude as in the present case [108]. State of the art experimental approaches for volumetric measurements in 3D unsteady flows are constantly evolving with new developments in both hardware and processing algorithms [116,117,[152][153][154][155]. Future technical developments providing increased spatial resolution should allow a more precise estimate of the velocity spatial gradients inherent to the viscous shear stress calculation, while preserving a measurement volume large enough to prescribe proper boundary conditions for the pressure gradient integration approach.

Euler-Bernoulli beam theory
As a comparison, Euler-Bernoulli beam theory is applied in exp. 1 to derive the transversal force distribution on the synthetic fin. This framework requires specific assumptions: that the deflections remain small (θ 0 = 10˚at the peduncle), the cross-sections do not deform and remain perpendicular to the deforming axis, and the deformations are confined in the proximo-distal direction. Under these assumptions, the hydrofoil obeys a linearized Euler-Bernoulli equation of motion, where the force per unit length normal to the surface (f n ) is given by [43]: Despite the fact that a quadratic function offers a good fit to the deflection profiles as described above, for a comparison with Euler-Bernoulli theory with stress-profiles another function h(x, t) has to be used. This is because a quadratic function would lack a non-zero fourth order spatial derivative and hence would not show a force profile. Therefore the timedependent deflections of the hydrofoil were once more fitted based on the raw images, but this time, with a function inspired by the general solution for a freely vibrating beam, clamped at one end, which is given by: The first term on the right hand side of Eq 12 is the inertial force, the second term, the elastic restoring force, the third term, the linear fluid damping and the fourth term, the internal viscoelastic damping, with m(x), the foil density per unit length, E, Young's modulus, I(x), the second moment of area, α and η, the fluid and internal viscous damping coefficients per unit length, respectively. The foil thickness d, its width B(x) and its density ρ f are used to evaluate m (x) and I(x): The boundary condition at x = 0, where the pitching motion θ(t) is imposed, can be expressed as: The boundary conditions at the free tip are not explicitly taken into account, since rather than solving Eq 12 for h(x, t), we use the known deflection to calculate the time-dependent load on the fin.
The fluid damping coefficient is determined based on the decay exponent γ obtained from the impulse response tests (Eq 1): This yields values of 0.009 kg m -1 s -1 in air and 0.03 kg m -1 s -1 in water, with a ratio of roughly 3.3 between both. An approximate ratio of 20 was found between α water and α air in a previous study, using a similar setup and a polysiloxane foil (with comparable thickness but larger area and larger Young's modulus) [43]. Based on a modal analysis of the beam deflection, a mathematical relation can be established between η, α, γ and f 0 (see equations B3 and B4 in [43]). Inserting our measured values of f 0 and γ (in air) into that equation, we estimated an internal damping coefficient (η) for our PDMS fin of the order of 10 -10 Nm 2 s. The internal damping contribution in Eq 12 is thus expected to be negligible compared to the other terms.

Control volume analysis
We used this approach as a complementary calculation to derive the time-dependent total force and verify the coherence of the PTV-based stress data. The control volume was defined as the whole 3D measurement domain, except the last two grid planes on each outer boundary walls.

Uncertainty analysis
A central point of our work is the evaluation of uncertainties on the hydrodynamic stresses, to ensure that stress maps can be measured accurately on the surface of a small synthetic fin with biologically relevant kinematics. To determine the positional uncertainty arising from the triplet reconstruction scheme, we consider the coordinates of a particle match, where the physical coordinates X i of the dewarped particles images are averaged over the different cameras N cam . From the uncertainty of each X i given by several parameters but dominated by the least square errors (LSE) from the calibration and the 2D Gaussian fit during the particle identification step, we can then obtain an uncertainty in position σ x using error propagation given by: The uncertainty in the z-direction (σ z ) is computed using the reference plane distance from the face plate of the cameras mount (*475 mm) and the distance between the apertures for a half-angle of about 6.5˚, leading to σ z ' 32μm. The positional uncertainty propagates to the velocity components: The time delay δt between two laser pulses is considered as exact, resulting in s u x ¼ s u y ' 0:002 m=s and s u z ' 0:018 m=s. During the velocity interpolation onto a regular 3D grid, appropriate filters were applied (see section 2.1). Therefore, these uncertainties are a conservative overestimation.
Noise propagation from the velocity field to the material acceleration and to the integrated pressure field has been the object of many studies [69,96,100,103]. Based on statistical arguments, a normal probability distribution can be determined for the pressure uncertainty, assuming a Lagrangian material acceleration, a normal probability distribution for the velocity uncertainty N ðm u ¼ 0; s 2 u Þ and integration paths covering n nodes [82]: Accordingly, the pressure error in the case of a zero-expectation value (μ u = 0) should not exceed the standard deviation: The Eulerian method would produce less clear error distributions due to the nonlinear advective term of the material acceleration in Eq 7, making it hard to define the uncertainty without a priori knowledge of the flow velocity field [82]. The Lagrangian approach is preferable to the Eulerian one as far as it allows a more direct treatment of the measurement parameters in the context of noise reduction. It was also demonstrated in other studies that the Eulerian approach suffers more from the velocity error propagation than the Lagrangian one, especially for flow fields characterized by strong advection velocity [97].
The main sources of uncertainty on the Lagrangian material acceleration are the truncation error associated to the numerical discretization, the precision error propagated from the velocity data and the pseudo-tracking scheme inaccuracy [97,103,158]. The discrepancy between the true and estimated particle positionsx p andx � p in the Lagrangian pseudo-tracking approach results in a discrepancy between the real and approximated velocitiesũ p andũ � p which adds up to the measurement uncertainty: Eqs 21 and 22 can be simplified by using the following approximation: A forward difference scheme is used to evaluate the material acceleration based onũ p ðtÞ andũ � p ðt þ DtÞ, yielding a total uncertainty on the material acceleration (with a i = Du i /Dt): It contains the error propagated from the velocity field and a second contribution from the path line reconstruction inaccuracy, whereas the truncation error from the interpolation of the velocity data, associated to the spatial resolution, was considered negligible compared to the other terms.
The uncertainty on the pressure is proportional to the uncertainty on the material acceleration, the spatial resolution Δx i and the number of nodes n crossed along the integration path. Moreover, the algorithm involves a median polling among a collection of eight paths. Since it is impossible to predict which direction will correspond to the selected median, averaging the uncertainty in all three directions is a reasonable approach. Besides, given a variable X with a statistical sample of size N and a median valueX, a mathematical relation between the variances is obtained [159]: The path length n � Δx i typically covers half of the domain size, which is of the order of 20 mm in the x-y plane, and 8 mm in the y-z plane. The spatial dependence of pressure errors is analyzed in section 3.4.
The discrete formulation of the viscous shear stress (Eq 5) is associated to an uncertainty: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Furthermore, the uncertainties are reduced by an additional factor of ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi N time � N spatial p , where N time is the number of time frames involved in the period-averaging step and N spatial is the number of grid points used in the spatial average. For the instantaneous stress curves of exp. 1, spatial averaging is performed over the most distal portion (15%) of the fin sides, covering the last four rows of points (N spatial = 4), whereas the stress curves extracted on the tip surface are obtained by averaging over the three rows of points covering that surface (N spatial = 3). In the case of the period-averaged distributions, N spatial = 4, as each stress curve is averaged over the left and right sides and over the dorsal and ventral portions of the foil.
Finally, the pressure and velocity uncertainties are used to estimate the error on the total force obtained from the control volume analysis, based on a discrete version of Eq 17, yielding the error bars of Fig 13.

Temporal and spatial resolutions
The time step between subsequent particles images (δt) should be fixed so that the particles do not move over more than *10 pixels during that time in the 2D captures (recommendation for the V3V setup from TSI Inc.). This ensures that the particles remain traceable by the algorithm, while limiting the error propagation to the velocity field. In our case, δt was set to 2.5 ms to respect that constraint. A similar twofold effect is associated with Δt and Δx, which must remain small enough to permit the Taylor expansions inherent to the calculation of the viscous shear stress and Lagrangian acceleration, but not too small to cause unreasonably large errors propagated from the velocity field. Additionally, the reliability of the pseudo-tracking approach underlying the Lagrangian acceleration can become an issue in the case of strongly curved trajectories involved in highly rotational fluid [69,103]. The simple criterion Δt ⩾ Δx/U ref was suggested [82]. The velocity magnitude of the characteristic flow structures is of the order of 0.075 m/s, implying a ratio Δx/U ref of about 0.01 s for exp. 1 (it is even smaller for exp. 2 where the vortical motion is stronger). Therefore, the condition stated above is met with a temporal resolution of Δt = 0.0125 s. Another criterion for vortex dominated flows was proposed, stating that a vortex should not execute more than half a turn during Δt, with a recommended a ratio of acquisition frequency (f acq = 1/Δt) to turnover frequency of at least 10 [97]. In our study, the acquisition frequency is 80 s -1 and the turnover frequency can be defined asũ ave =a, yielding values of 1.4 s -1 and 2.07 s -1 for exp. 1 and 2, respectively, which shows that the criterion is fulfilled.
For the queen2 algorithm, pressure errors in the range of 5-10% were reported for Δx < 0.0625a, where a is the characteristic dimension [77]. The flapping foil in our experiments covers a displacement of 11-14 mm, which translates into a ratio of Δx/a between 0.054 and 0.068, indicating that the grid spacing is small enough to produce reasonably small pressure errors. Lastly, a recommendation was made to ensure the validity of the null-pressure boundary condition inherent to the queen2 algorithm, based on the condition H/a ⩾ 2 where H is half of the domain size [77]. Substituting the tip excursion amplitudes (11-14 mm) for the characteristic dimension a and replacing H by the average half domain size in the x and y directions, we conclude that the size of the measurement domain lies right above that limit.

Hydrodynamic stress distributions
3.1.1. Experiment 1. The parametrized fin deflections are presented in Fig 4 (right panel), at eight selected time points for which the stress distributions are presented. The most proximal portion of the fin is excluded from the 3D particle position field (Fig 4, left panel), and in the following graphs, the zero of the x axis is relocated at the initial point on the fin where the flow quantities are extracted.
The velocity field at t 0 is shown in Fig 5, where 2D slices from the measurement volume are presented: one x − y plane (in the middle of the domain) and two x − z planes located on both sides of the fin. In all planes, the shedding of counterrotating vortices can be observed, corresponding to the 2D projections of vortex rings. In Fig 6, the velocity field at t 1/4 is represented by a color map. The velocity magnitude is depicted in the x − y midplane, where high values are associated to two vortices, one previously shed in the downstream wake and a second still in formation, attached to the foil. The vertical (u y ) and transversal (u z ) velocity components in a single x − z plane are also shown in Fig 6. A pair of subsequent time points (t 1/4 and t 3/8 ) is analyzed to ensure that the stress profiles on the fin surface present a realistic time-evolution (no unphysical abrupt transitions). The left side of the foil at t 1/4 is subjected to positive pressure, yielding a normal stress directed towards its surface. The maximum value is of the order of 10 Pa, located distally on the fin, although not at the very end, and centrally along the dorso-ventral axis. The tip corners experience negative pressure with values around -4 Pa. A similar trend is observed for the left side in the following time frame (t 3/8 ), where the pressure maximum is reduced to a value around 8 Pa. At t 1/4 , the right side of the fin is subjected to negative pressure on its most distal portion, of the order of -7 Pa, an effect that concentrates more centrally along the dorso-ventral axis at t 3/8 , with a minimum decreasing to approximately -15 Pa, along with the appearance of two small positive pressure spots of about 5 Pa close to the dorsal and ventral corners.
The relative contribution of the viscous stress within the total hydrodynamic stress, in comparison with the pressure, can be estimated based on the distributions of the six components τ ij from Eq 4, namely, the viscous stress acting in the direction i on a surface perpendicular to j. These six components are presented in the x-y midplane for t 1/4 in Fig 8. As expected for Reynolds numbers corresponding to the inertial flow regime, these distributions are about two orders of magnitude lower as compared to the pressure. Local variations on the virtual fin surfaces are visible, corresponding to the viscous shear profiles within the boundary layer. As explained in section 2.4, the cartesian components (τ xx , . . ., τ xz ) are subsequently projected on the surface unit vectors and decomposed into proper components in the reference frame of the fin (proximo-distal, dorso-ventral, left-right). Apart from these noticeable variations in the vicinity of the solid surface (especially for τ xx , τ yy and τ xy ), special flow signatures are observed within the vortex that was identified from its negative pressure core in Fig 7, including paired regions of opposite-signed viscous stress, typical of the rotational nature of the vortex. For an incompressible fluid such as water, the divergence-free condition (r �ũ ¼ 0) implies that the sum of τ xx , τ yy and τ zz should be equal to zero. This condition was verified for the instantaneous flow field presented in Fig 8, where the divergence of the velocity vector field was found to be within ±15 s -1 , with an average over the whole volume (in absolute values) of 2.4 s -1 . The uncertainty on the divergence ofũ can be evaluated using error propagation from the velocity components, yielding a value of 17.2 s -1 (with centered finite differences). Therefore, the measured flow field is divergence-free within experimental uncertainties.
The examination of instantaneous stress maps at mirroring time points (equivalent motion instants such as t 0 and t 1/2 or t 1/4 and t 3/4 ) is used to corroborate the consistency of the symmetrical distributions. Spatial symmetry between the left/right and ventral/dorsal sides of the hydrofoil is also verified to validate the stress maps. The stress patterns are drawn more explicitly by plotting the individual components along the relevant axes. The curves shown in  represent the normal stress component (at t 1/4 and t 3/4 ) plotted against the dorso-ventral axis, averaged over the most distal portion (last 15%) of the membrane. By comparing both panels of Fig 9, we notice that as the fin is passing through mirroring sites, the normal stress on the left and right sides are mirrored. The absolute value of the pressure is maximal at the center of the foil and decreases towards the dorsal and ventral edges. The side which is leading the motion (the left side at t 1/4 and the right side at t 3/4 ) experiences a negative normal stress, which means that the stress is directed towards the surface. On the the trailing side, the positive normal stress points in the direction of the outward normal, as the fluid tends to hold the foil back. Small discrepancies between the mirroring time points can be attributed to experimental errors perturbing the symmetry of the foil motion, such as an imperfect positioning relative to the incoming flow.
In Fig 10, we consider the viscous shear stresses on the tip surface, at mirroring time points t 0 and t 1/2 . At t 0 , the fin travels towards its right side (positive y values) and the tip surface experiences shear stress pointing in the left direction (S lr < 0), with a maximum absolute value in the center. The opposite trend is observed at t 1/2 . As for the dorso-ventral shear stress, it is positive on the ventral portion of the tip, and negative on the dorsal part, implying that the stress is directed towards the center as the fluid tends to contract the fin tip surface. Naturally, this effect occurs for both t 0 and t 1/2 since the dorso-ventral axis is symmetric between the mirroring time frames.
Finally, the period-averaged absolute value of the normal stress component (averaged over a full cycle) is illustrated in Fig 11. Three curves are shown, each one corresponding to a selected position along the proximo-distal axis of the fin: at 25%, 50% and 100% of the total length (0.25L, 0.5L and L).

Experiment 2.
An interesting application of the hydrodynamic stress mapping method lies in the investigation of particular stress signatures or thresholds along the main axes of the hydrofoil, where the effects of varying specific flow parameters can be analyzed. The period-averaged absolute values of the normal stress are compared between exp. 1 ( Fig  11) and exp. 2 (Fig 12) to deduce how the hydrodynamic load is changing over the fin surface depending on the flow conditions (without or without an external current). These distributions correspond to averages over complete cycles. Experimental imperfections in the foil geometry and its positioning inside the flow chamber can result in a slight asymmetry of the fluid force distributions on the dorsal and ventral sides. Averaging the stress curves over both sides allows to reduce noise and experimental errors. For the distributions prior to dorso-ventral averaging, the error bars are not shown for sake of clarity, as they would be larger than for the symmetrical curves by a factor of ffi ffi ffi 2 p . The magnitude of the normal stress at 50% of the fin length (0.5L) is reduced as a consequence of applying an external flow (exp. 2). Besides, at the distal location on the foil (L), the curve goes from a single peak for exp. 1 to a curve with two maxima for exp. 2, with similar maximum magnitudes. Therefore, the presence of an external advecting flow, combined with slightly larger flapping amplitude, shifts the stress curves to lower values at specific positions along the fin (about half-way through the proximo-distal axis) and modifies the peaks profile on the most distal surface. These results support the practicality of the stress mapping method as a powerful tool for analyzing the distribution of fluid forces on submerged and deforming objects.

Control volume analysis
A control volume analysis was used to validate the consistency of the instantaneous force data. The total forces acting on the fin in the x and y directions (F x and F y ) were calculated from Eq 17, for the eight selected time points of exp. 1. The results are shown in Fig 13. Fair agreement is found between the control volume forces and the values obtained by integrating the stress distributions on the fin surfaces within the uncertainty of the pressure determination corresponding to an uncertainty in the force of σ F = 0.4 mN for the PTV-integrated results and σ F = 0.2 mN for the control volume method. Both approaches describe a net force in the negative x direction with a maximal magnitude of *0.75 mN at mirroring instants t 1/4 and t 3/ 4 (where the hydrofoil tip excursion is maximal), implying that the foil experiences a net thrust as it is propelled upstream. Both methods also capture a net force in the positive y direction in the first half of the period (instants t 1/8 to t 3/8 ) and in the negative y direction in the second half of the period (instants t 5/8 to t 7/8 ), with a maximal magnitude of *1.75 mN. We conclude that the time-evolution of the PTV-based stress distributions is consistent with the fluid forces obtained from a control volume analysis. A small offset is seen in F y at t 3/8 and t 7/8 , which can be attributed to experimental errors in the velocity field. To investigate the source of this offset, we performed again the control volume calculation, using a virtual volume slightly closer to the fin in the x and y directions (not in the z direction though to ensure that the whole fin remains included in the control volume). The results are identical within the uncertainties, indicating that the shift originates from the PTV-based stress integrated over the fin surfaces, and not from the control volume. Globally, the magnitude and shape of the time-resolved total forces are in good agreement between both approaches.

Euler-Bernoulli force distributions
The motion instants t 0 , t 1/8 , t 1/4 and t 3/8 from exp. 1 (covering half a period) were selected to evaluate transversal forces in the Euler-Bernoulli framework. The parametrized deflection of the foil (Eq 13) was used to compute the appropriate derivatives in Eq 12, and to provide an estimate of the transverse force per unit length (pointing upwards) along the proximo-distal axis of the fin. The results are presented in Fig 14, and compared with the force distributions obtained by integrating the PTV-based normal stress distributions over both sides of the hydrofoil. The main contribution to the total transverse load is the elastic restoring force, which is of the order of 0.1 N/m, followed by the inertial force, of the order of 0.01 N/m. The linear fluid damping and the structural damping contributions are negligible in comparison, of the order of 0.001 N/m. The curves from both approaches are consistent in terms of order of magnitude and proximo-distal trends. The force drop at the tip of the fin for t 1/8 and t 1/4 , captured by the PTV-based stresses, is not well reproduced by the Euler-Bernoulli model. We attribute this small discrepancy to the larger tip deflections at these specific time points, where the Euler-Bernoulli formalism is less applicable by definition. Moreover, the quality of the fit (Eq 13) is slightly decreased at the tip of the fin for these time points, where the deflection is overestimated by the fitted function. As expected, the Euler-Bernoulli calculation offers a more approximate solution to the force maps, where the hydrodynamic forces obtained directly from the 3D flow field provide more detailed local stress variations. In brief, the force distributions obtained from the Euler-Bernoulli deflection analysis compare reasonably well with the PTV-based normal stresses, given the underlying assumptions.

Uncertainties
From Eq 20, it is possible to provide a raw estimate of the pressure uncertainty, using the parameters specific to our experiments. With an approximate integration path length (covering half the fluid domain) of 30 nodes in the x and y directions and 10 nodes in the z direction, we get an approximate pressure error of 1 Pa (after averaging over the x, y, z directions).
A more detailed analysis is rendered possible by the use of Eq 27. In the case of exp. 1, the error on the material acceleration arising from the Lagrangian path reconstruction inaccuracy (the second term in Eq 25 which carries the spatial dependence) was calculated everywhere in the 3D domain. Spatial average over the whole volume and temporal average over a complete period (using the 8 time frames t 0,. . .,7/8 ) yielded the following values:

PLOS ONE
The larger value obtained for the z component contributes to larger uncertainties accumulating in that direction. In comparison, the first term of Eq 25, which contains the material acceleration uncertainty propagated from the PTV errors, produces values at least one order of magnitude larger: We thus conclude that the error propagated from the velocity field constitutes the main source of error in the pressure calculation, whereas the contribution from the inaccuracy of the Lagrangian path reconstruction is almost negligible. The resulting pressure uncertainty (σ p ), averaged over the 3D domain and a complete motion cycle, is equal to 4.03 Pa. The spatial dependence of σ p is shown in Fig 15 (for t 0 , exp. 1). The fluid vorticity (õ) and its streamwise component (in absolute value, ω x ) are illustrated in Fig 16. The pressure uncertainty is found to be spatially correlated to the vorticity, in particular, the error becomes larger vis-a-vis the tip corners, where high vorticity around thex axis is found. This correlation is expected since the spatial dependence of σ p is carried by the term associated to errors in the pseudo-tracking method. Thus, it is no surprise that the regions where the particles follow paths with strong curvature (high vorticity) are the regions where the Lagrangian path reconstruction, which assumes a linear displacement between subsequent velocity fields, performs more poorly.
Including the reduction factor due to spatial and temporal averaging, the final uncertainties on the normal stress in exp. 1 are σ p ' 2 Pa for the instantaneous distributions (Fig 9) and 0.7 Pa for the period-average (Fig 11). The uncertainties on the viscous shear stresses (acting on the tip, Fig 10) are σ τyx ' 1.6 mPa for the left-right component and σ τzx ' 10 mPa for the dorso-ventral component. However, the dorso-ventral symmetry and the temporal symmetry in the shear stress curves suggest on the contrary that these stress distributions are well resolved within an uncertainty range that does not exceed 3 mPa. We interpret this as an indication that the velocity errors stated in section 2.7 are slightly overestimated, as they do not include the appropriate error reduction associated to the filtering and smoothing steps. For that reason, the error bars were exceptionally omitted on that specific graph (Fig 10, lower  panel). For exp. 2, the error bars in Fig 12 correspond to σ p ' 0.35 Pa (in that case, the periodaverage was performed over 32 frames).
Finally, the error bars of Fig 14 were calculated based on σ p , the integration path length (the average fin width) and the number of grid points along that path (25 on each side of the fin).

Discussion and outlook
In the present study of a prototypical fin-like foil, we performed three-dimensional hydrodynamic stress calculations based on a pressure gradient multi-path integration scheme, which we applied to a time series of volumetric PTV measurements of an unsteady vortical flow. This method was shown to produce well resolved stress curves along the biologically relevant axes of the synthetic fin, allowing to gain precise information about the interaction between the deforming surface and its direct fluid environment. We explored in details the error propagation for both the normal and tangential stresses. We validated the results using a control volume analysis and compared the force trends with values determined from the Euler-Bernoulli beam theory.
The non-intrusive approach of hydrodynamic stress mapping constitutes a promising path to deepen our understanding of the relation between the morphometric and elastic properties of hydrofoils and the fluid forces which they experience. We foresee numerous applications to this experimental methodology, in domains such as the engineering of underwater propulsive foils or the investigation of hydrodynamic forces in aquatic animal locomotion. Propulsive foils involved in marine robotic systems would greatly benefit from the employment of small scaled models for which the distribution of hydrodynamic stresses and their time evolution could be described everywhere on the surface, without the need for intrusive mechanical transducers, thus guiding the structural design with precise knowledge of the local mechanical loads. In a parallel train of thought, mapping the stress patterns on the surface of biomimetic synthetic fins can provide new insights about the mechanical constraints under which real fins have evolved to attain their exact morphology.