Connections between propulsive efficiency and wake structure via modal decomposition

Abstract We present experiments on oscillating hydrofoils undergoing combined heaving and pitching motions, paying particular attention to connections between propulsive efficiency and coherent wake features extracted using modal analysis. Time-averaged forces and particle image velocimetry measurements of the flow field downstream of the foil are presented for a Reynolds number of $Re=11\times 10^3$ and Strouhal numbers in the range $St=0.16\unicode{x2013}0.35$. These conditions produce 2S and 2P wake patterns, as well as a near-momentumless wake structure. A triple decomposition using the optimized dynamic mode decomposition method is employed to identify dominant modal components (or coherent structures) in the wake. These structures can be connected to wake instabilities predicted using spatial stability analyses. Examining the modal components of the wake provides insightful explanations into the transition from drag to thrust production, and conditions that lead to peak propulsive efficiency. In particular, we find modes that correspond to the primary vortex development in the wakes. Other modal components capture elements of bluff body shedding at Strouhal numbers below the optimum for peak propulsive efficiency and characteristics of separation for Strouhal numbers higher than the optimum.


Introduction
The performance characteristics of swimming and flying animals have long motivated the design of autonomous swimmers with similar kinematics for propulsion and navigation.The energetics of fish propulsion and, in particular, their ability to harvest energy in schools or unsteady flows environments (Beal et al., 2006) have also spurred the research and development of novel energy harvesters (McKinney and DeLaurier, 1981;Jones et al., 1997;Bryant et al., 2012).Broadly, these animals rely on oscillating or undulating foils for propulsion.
Characterizing the relationship between foil performance and the downstream wake structure has been a long-standing research area.Pioneering work in this area by Triantafyllou et al. (1993) analyzed the swimming performance of several different fish and showed that fish operate within a narrow Strouhal number range, 0.2 < St < 0.35, that leads to high propulsive efficiency.The Strouhal number is a kinematic parameter defined as St = f A/U ∞ , where A is a characteristic length describing the width of the wake, f is the frequency of oscillation, and U ∞ is the swimming speed.Laboratory experiments also showed high propulsive efficiencies from oscillating foils operating in this Strouhal number range.Complementary stability analyses showed that oscillating foils produce a jet-like wake that is convectively unstable when excited at frequencies corresponding to Strouhal numbers that yield peak propulsive efficiencies.The spatial stability analysis performed by Triantafyllou et al. (1993) used the jet profile of a reverse von Kármán vortex street, where two single (S) vortex cores were shed per half cycle of oscillation, i.e., a 2S wake pattern.Similar relationships between optimal Strouhal number ranges and instabilities of the jet wake were found by Lewin and Haj-Hariri (2003).
The relationship between wake instability and peak propulsive efficiency was studied further by Moored et al. (2012), who used Particle Image Velocity (PIV) to visualize wake structures produced by a batoid-inspired oscillating fin.These experiments showed multiple peak efficiencies, with some corresponding to a 2P wake pattern, where two pairs (P) of vortex cores are shed per half cycle, and others corresponding to a 2S wake.Moored et al. (2012) also pursued local stability analyses (i.e., relying on a parallel flow assumption) to show that the time-averaged velocity profile exhibits instabilities when excited at frequencies corresponding to peak propulsive efficiency.This phenomenon was termed 'wake resonance', suggesting that when the oscillating foil is tuned to optimally excite the wake, it produces the highest propulsive efficiencies.Moored et al. (2012) also examined the vorticity perturbations generated by the most unstable modes in an effort to better understand these instabilities.Follow-on work by Moored et al. (2014) showed that the 'resonant frequencies' correspond to optimal momentum entrainment, both into and out of the jet region.These findings suggest that there is no single wake structure that corresponds to peak propulsive efficiency (Smits, 2019).Further examples of 2P vortex patterns have also been found in wakes produced by eels and dolphins (Tytell and Lauder, 2004;Smits, 2019).
Recent studies have raised questions regarding the validity of wake resonance theory (Arbie et al., 2016), noting that although correlations exist between unstable frequencies and peak propulsive efficiency from experiments, it is more difficult to establish a causal link between the two.Arbie et al. (2016) considered the stability characteristics of momentumless wakes and noted that these wakes may be stable (even if the thrust producing jet-like component is unstable).It has also been suggested that, while wake structure provides insight into propulsive efficiency, it cannot provide a complete explanation (Zhang, 2017;Taylor, 2018).Other studies (Eloy, 2012) suggest that the development of the wake structure is a result of, and not a cause of, high propulsive efficiency.The reverse von Kármán wake from an oscillating foil can be shed in more patterns than just the 2S and 2P wake (see e.g., Lentink et al., 2008;Schnipper et al., 2009;Andersen et al., 2017), and each wake structure has an indirect relationship with propulsive efficiency.For example, Mackowski and Williamson (2015) found from experiments that wake patterns for a pitching foil, where self-interactions of the vortices occur, do not reflect a net force.Thus, observing the wake structure alone is likely not enough to completely determine the propulsive efficiency (Floryan et al., 2020).
As an alternative approach, we make use of dynamic mode decomposition (DMD) to identify coherent features (or patterns) in the wakes produced by oscillating foils and link these features to propulsive performance.DMD, as first introduced by Schmid (2010), is a technique used to approximate the dynamics of a nonlinear system via the identification of a linear operator that evolves the system to the next state.DMD can also be thought of as a data-driven stability analysis because each DMD mode is associated with a specific frequency and growth or decay rate, which provides interpretable physical insight into the spatial structures and their dynamic evolution.In contrast to the stability analysis approach, DMD does not require the mean flow to be locally or globally parallel.In the fields of marine propulsion and energy harvesting, DMD and similar modal decomposition techniques have been used on propeller and turbine wakes to characterize wake instabilities, loading conditions, and efficiency (Magionesi et al., 2018;Sarmast et al., 2014;Strom et al., 2022;Araya et al., 2017).
In this study, we use the triple decomposition to identify coherent flow features that contribute to drag and thrust production from PIV measurements for the wake past an oscillating foil.We examine both the vorticity and the Reynolds stresses associated with optimized DMD modes across Strouhal numbers ranging from St = 0.16 to St = 0.35, and compare the results to time-averaged forces and propulsive efficiencies.We find that the use of DMD enables us to link wake structure with wake stability and propulsive efficiency for oscillating foils.
Our study builds upon the foil configuration described by Floryan et al. (2020).We focus on rigid foils and present results primarily in the context of swimming performance.While the use of rigid foils represents a simplification of the fluidstructure interactions pertinent to flapping foil propulsion in nature, our results may provide additional insight into coherent flow features that contribute to drag and thrust.Moreover, the observations presented herein are also relevant to the design of autonomous swimming vehicles (Buren et al., 2020) and energy harvesting systems (McKinney and DeLaurier, 1981) that use rigid oscillating foils.

Heave and pitch foil parameters
We consider an oscillating NACA-0012 hydrofoil with chord c and span s as illustrated in figure 1.The distance b = 0.25c denotes the point of rotation from the leading edge.The imposed pitching and heaving kinematics of the foil at the point of rotation are described by θ(t) = θ 0 sin(2π f t + ϕ p ) and h(t) = h 0 sin(2π f t) respectively, where θ 0 and h 0 are the respective pitching and heaving amplitudes, f is the frequency of oscillation and ϕ p is the phase difference between heaving and pitching oscillations.For the experiments described below, the phase difference was set to π/2 (90 • ).
The oscillation frequency can be expressed in dimensionless terms as the Strouhal number: Note that the Strouhal number can also be interpreted as the wake width (∼ A TE ) normalized by the wavelength (∼ U ∞ / f ).Another important parameter relative to swimming performance is the angle of attack α, which is different from the pitch angle due to the influence of the heave velocity, ḣ(t), and can be expressed as: (2) In the case where the phase angle ϕ p between pitch and heave is 90 • , the maximum angle of attack α 0 can be estimated as: where ω = 2π f is the radian frequency.The time averaged power and thrust coefficients are described using the following equations: where ρ is the density of the fluid, τ is the time-averaged thrust, and ℘ = dt is the time-averaged power input to the fluid.F y and T z are the force and torque associated with the heaving and pitching motions.Propulsive efficiency can then be calculated as: The propulsive efficiency and thrust coefficient are both used to characterize swimming performance.

Experimental Methods
Experiments were carried out in a large-scale free-surface water channel with a glasswalled test section of dimensions 7.6 m × 0.6 m × 0.9 m.A NACA-0012 hydrofoil with a chord of c = 10 cm and a span of s = 32 cm was used.The foil was 3D printed from polylactic acid (PLA) filament, sanded, and reinforced with two internal aluminum rods throughout the span.The channel was run at a freestream speed of U ∞ = 0.1 m/s with a turbulence intensity of 1%.Foil motions were produced using a closed-loop control system driven by two NEMA 23 integrated stepper motors, as shown in figure 2. One of the motors was connected to a linear rail to generate the heaving motion.The other motor was mounted to the beam to create the pitching motion.Translational and angular positions of the foil were calculated from the precision micro-steps (3200 steps per revolution) of the stepper motors.Measurements of hydrodynamic forces were made concurrently with an ATI Gamma (SI-32-2.5)six-degree-of-freedom force transducer with a minimum force and torque resolution of 1/80 N and 1/2000 Nm respectively.Data were acquired at a rate of 5 kHz and filtered using a zero-phase second-order Butterworth filter.Instantaneous forces were measured over oscillation frequencies f = 0.4 Hz−1.3 Hz.For the kinematic parameters used in the present study (Table 1), this yielded a Strouhal number range of St = 0.16 − 0.59.Time-averaged thrust coefficients (C T ), power coefficients (C P ), propulsive efficiencies (η) were computed by combining the force measurements with the foil kinematic data.The standard Klein-McClintock procedure was used to estimate uncertainties for these parameters based on known instrument resolution and measurement uncertainty for the time-averaged thrust and power input (Taylor, 1997).Force measurements were collected and averaged for at least 13 periods for each run.Freestream velocity measurements were recorded using a Laser Doppler Velocimetry system (MSE miniLDV).Relative uncertainties in the velocity measurements were 0.11%.2D-2C PIV measurements were carried out in the near-wake region of the oscillating foil to provide further insight into the transition from drag-producing and freeswimming conditions to thrust-producing conditions.The PIV system comprised a 5 W 532 nm solid-state laser and a Phantom VEO 410-L high speed camera with 1280 × 800 pixel resolution.Images of the flow field were captured in a field-of-view of 27 cm × 17 cm at a rate of 65 Hz for approximately 56 s.Standard analysis routines in DaVis (LaVision) were used to generate velocity measurements.Two 64 pixel × 64 pixel and four 24 pixel × 24 pixel convolution windows were used with a 50%

Periodic Structures via the Triple Decomposition
For periodic flows under natural or forced conditions, large-scale coherent motions are often present in addition to turbulence.In such situations, the Reynolds decomposition, which assumes that turbulence is the only source of fluctuations in the flow, may not be accurate and can lead to overestimation of the stochastic part of the flow.Alternatively, the full flow field u can be expressed as a triple decomposition as introduced by Hussain and Reynolds (1970): where ū is the mean (time-averaged) flow field, ũ is the periodic flow field, characterized by the phase parameter ϕ(t), and u ′ represents the turbulent fluctuations that are incoherent.In our case, where the flow is driven by a known forcing frequency, the periodic component can be directly obtained through phase-averaging as: where the specific phase position ϕ 0 = ϕ(t 0 + nτ) is defined for an initial time t 0 and a period τ.However, for natural flows, or those driven by a series of forcing frequencies, equation 7 may be difficult to compute.Instead, phase averaging can be attained through modal decomposition methods.In these cases, the unsteady Figure 3: Hydrofoil and PIV experimental setup with field-of-view dimensions.
coherent component ũ may be represented by a linear combination of modes obtained via techniques such as Proper Orthogonal Decomposition (POD) or Dynamic Mode Decomposition (DMD), and the triple decomposition can be modified as described by Baj et al. (2015): Here, r is the number of modes retained from the modal decomposition.In this study, we use Dynamic Mode Decomposition to obtain the periodic component of the flow field in equation 8 and compare this modal representation to the wake characteristics obtained via direct phase-averaging, as shown in equation 7.

Dynamic Mode Decomposition
DMD was originally presented to the fluid mechanics community by Schmid (2010) as a method to decompose unsteady or turbulent flow fields into coherent structures.DMD can be interpreted as an eigendecomposition of a least squares best-fit linear operator A that advances the past snapshots of data (i.e., 2D-2C velocity fields obtained from PIV) forward in time towards future snapshots.For a set of m snapshots (x 1 , x 2 , • • • , x m ), the dataset is arranged into two separate matrices of the following forms: and Generally, A can be approximated by taking the Moore-Penrose pseudo inverse of X ′ (A = X ′ X −1 ) via the singular value decomposition (SVD): X = UΣV * .However, for high-dimensional datasets such as measurements from PIV snapshots, the computation of A is expensive and can be inaccurate due to noisy entries or outliers in the data.To limit computational expense and the impact of noise, we consider the projection of A onto a reduced-rank representation of the data obtained via the SVD, X ≈ U r Σ r V * r , where the subscript denotes a truncation to rank-r: Note that this is analogous to projecting A onto the leading POD modes.Subsequently, the eigenvalues of matrix A can be estimated by taking the eigendecomposition of Ã: where the matrix W contains the eigenvectors and Λ contains the corresponding eigenvalues (λ 1 , λ 2 , • • • ).Dynamic modes for the flow field can be computed using the exact DMD algorithm from Tu et al. (2014) as: Here the matrix Φ contains individual dynamic modes (ϕ 1 , ϕ 2 , • • • ).Using this lowdimensional approximation to the dynamics represented by A, we can reconstruct the data using a linear approximation of the dynamic modes for all future times: The frequency ω k is defined based on the associated eigenvalue λ k as ω k = ln(λ k )/∆t, where ∆t is the time interval between individual snapshots.The vector b contains entries of the initial amplitudes b k defined as b = Φ −1 x 1 .Recall that the dynamic modes ϕ k are the columns of the matrix Φ.The modes in this temporal DMD formulation represent the absolute stability of the flow field, with the eigenvalues providing insight into oscillation frequencies and growth or decay rates.Given a complex eigenvalue , the frequency can be estimated as ∆t .Thus, a dynamic mode ϕ k grows over time if the corresponding eigenvalue has amplitude |λ k | > 1 and decays if |λ k | < 1.
Spatial growth rates can also be approximated by the DMD algorithm by reorganizing the series of snapshots as increasing in space.In this case, modes in spatial DMD represent convective stability, where the eigenvalues would represent spatial frequencies or wavenumbers.In practice, spatial DMD is more prone to noise due to the sparsity of spatial data from PIV measurements, since the number of timesnapshots available is typically much higher than the spatial locations (Schmid, 2010).For instance, the PIV data obtained in this study span 3600 snapshots in time but only 77 streamwise locations for each run.As a result, we focus on using the temporal DMD approach to take advantage of our time-resolved snapshot data.

Optimized DMD
For the periodic flows that are studied in this paper, it is expected that the dynamic modes should neither grow or decay in time.That is, the eigenvalues of DMD are fixed directly onto the unit circle |λ k | = 1.However, the presence of even weak noise in periodic flows is known to yield damped eigenvalues, |λ k | < 1 that result in decay of the corresponding modes (ϕ k ) over time (Bagheri, 2014).Thus, several recent methods have been proposed for debiasing the DMD algorithm in the presence of noisy data, such as measurements from PIV (Dawson et al., 2016;Hemati et al., 2017;Askham and Kutz, 2018).For this purpose we use the optimized-DMD (opt-DMD) algorithm presented by Askham and Kutz (2018).This is a variant of DMD that uses a variable projection method to approximate the linear operator A with reduced noise.Specifically, opt-DMD solves the non-linear least-squares problem: where the eigenvalues in Λ act as iterative parameters and determine the values of the eigenfunctions, Φ, and amplitude coefficients, b.The algorithm also has the advantage of projecting the modes onto the full dataset rather than a subset from an r-rank truncation.Additional constraints can be applied to the modes, such as restricting eigenvalues to stay on the unit circle.For simplicity, we do not use any constraints for the flow fields in this study, and find that opt-DMD naturally converges towards the expected eigenvalues in a periodic flow.The DMD modes of velocity are computed using both opt-DMD and exact DMD algorithm for comparison.We compute the vorticity directly from the opt-DMD modes of velocity to understand the differences between the modal representation and the full flow field.

Propulsive Efficiency
Measurements of the thrust coefficient and propulsive efficiency are presented in figure 4. The thrust and power coefficients increase monotonically with increasing Strouhal number.However, C P increases faster than C T at higher St values, leading to non-monotonic behavior in propulsive efficiency and reduced efficiency for St > 0.35.Peak propulsive efficiencies can be seen within the range of 0.23 < St < 0.35.With the provided uncertainty estimates, it can be argued that either of the three Strouhal numbers St = 0.23, 0.29, and 0.35 result in the highest propulsive efficiency.However, simplified scaling laws suggest that for angles of attack α where the flow remains attached to the foil, peak propulsive efficiency can be estimated as where St 0 is the Strouhal number that results in zero net thrust (Taylor, 2018).In our case, the thrust coefficient is near zero for St = 0.16.Using this value for St 0 gives a maximum efficiency of St max ≈ 0.28.This approximation is consistent with our case of St = 0.29 for which the measured η is highest.Note that the uncertainty in measured efficiency at St = 0.16 is large, which is due to the fact that the thrust and power coefficients are close to 0 for this condition.Nonetheless, this case is similar to a self-propelled swimming mode, where the thrust from the foil is approximately equal to its drag.Negative efficiencies are expected for lower Strouhal numbers, St < 0.16, for which the net thrust becomes negative due to the effects of viscous drag.
As the foil oscillates at lower frequencies, the drag contribution on the propulsor stays approximately constant while the thrust generated decreases (Floryan et al., 2017).The present measurements are broadly consistent with trends observed from previous studies (Quinn et al., 2015;Triantafyllou et al., 1993).

Mean and Phase-averaged Flow Fields
PIV measurements were used to study the downstream wake of the foil for a subset of the Strouhal numbers for which the force measurements were obtained.Phaseaveraged vorticity (Ω) and time-averaged streamwise velocity ( ū) fields are shown in in figure 5, where the foil is positioned to move in the positive y direction.An increase in magnitude of the jet wake (i.e., ū/U ∞ > 1) can be seen with increasing St.A near-momentumless wake is observed for St = 0.16, where the wake profiles are close to zero, again indicating a resemblance to self-propelled swimming (figure 5b).Under this condition, approximately five positive and six negative vortices are shed per cycle of oscillation, signifying a small wake asymmetry.The vortices deflect laterally away from the heave and pitch centerline, resulting in an increase of the wake width from x/c ≈ 0.5 − 1.5 that can be seen in figure 5b.
Increasing the Strouhal number to St = 0.23 produces a 2P wake (figure 5c), with mean profiles that are approximately trapezoidal-shaped, with two small but distinct peaks.This wake characteristic matches well with the 2P structure observed by Dewey et al. (2012), who concluded that the presence of two vortex pairs result in two separate peaks in the mean profile.The Strouhal number St = 0.29 (figure 5e), initially produces a 2P wake with mean profiles that are trapezoidal (0.22 ≤ x/c ≤ 1.5).These vortex pairs coalesce further downstream to a 2S wake, creating a mean profile resembling that produced by a single jet (1.5 < x/c ≤ 2.7).As the Strouhal number increases further to St = 0.35, a classical 2S wake is observed, suggesting a transition from 2P to 2S wake structure over the range 0.23 < St < 0.35.These observations are consistent with the results of Moored et al. (2012), who found optimal efficiencies in both wake types.

Modal Decomposition
We now compare results between the exact-DMD and opt-DMD algorithms.Eigenvalues obtained for the first 9 DMD modes for both methods are shown in figure 6.
The mode corresponding to λ = 0 for each Strouhal number represents the mean flow field, while the other eigenvalues, having a non-zero imaginary part, represent the time-periodic modes.Note that the eigenvalues with non-zero imaginary components are represented in complex conjugate pairs (i.e., with λ k = λ k,r ± λ k,i ), which together represent a single frequency component of the flow field.
As shown by Magionesi et al. (2018), the eigenvalue problem degenerates and results in harmonic solutions when the dominant components of the flow travel at a constant speed, and with constant shape.This is noticeable in both the opt-DMD and exact-DMD modes, whereby the eigenvalues have oscillation frequencies (|ω k,i |) that are multiples of the first harmonic.The opt-DMD eigenvalues for all Strouhal numbers lie directly on the unit circle, suggesting that the modes neither grow or decay overtime (|λ k | = 1).In contrast, the exact DMD eigenvalues tend to fall within the circle for increasing St, particularly for large values of |λ k,i |, signifying that the oscillatory modes obtained via exact DMD decay over time (|λ k | < 1).This is most evident for the largest Strouhal number case (St = 0.35) in figure 6d, which shows eigenvalues within the unit circle for |λ k,i | > 0.2.Similar trends were reported by Bagheri (2014), who observed that eigenvalue damping increases linearly with noise amplitude and quadratically with frequency.For the remainder of this paper, we refer to the flow field produced by the dynamic modes associated with a complex conjugate pair of eigenvalues as a single mode, and consider the first to fourth harmonics as modes one to four.
To further understand the effects of mode decay via eigenvalue damping, we computed the periodogram-based power spectral density of the original flow field and the reconstructed flow from opt-DMD and exact-DMD modes (see figure 7).Note that each of the peaks represents the flow field produced by a pair of dynamic modes (i.e., corresponding to a complex conjugate pair of eigenvalues).The normalized frequen-  cies for these peaks, St = f A TE /U ∞ , represent harmonics of the foil oscillation frequency.The power from the exact-DMD reconstruction are lower in power compared to the opt-DMD reconstruction and full flow field (figure 7c and d).This decay can be attributed to the dampened eigenvalues from the exact-DMD algorithm, observed in figure 6(c,d).As a result, the effect of reduced power worsens with increasing St.These results suggest that the opt-DMD algorithm outperforms the exact-DMD method in retaining the energy in all modes, as also shown by Strom et al. (2022).
For completeness, we also highlight the exponential decay of the amplitude peaks of the opt-DMD modes, as shown in the dotted line of figure.Interestingly, the exponential fits for St = 0.16 and St = 0.23 show some flattening at the higher frequencies, suggesting that the wake is not strongly locked to the actuation frequency for these cases.For example, in the case of St = 0.23 (figure 7b), the third opt-DMD and exact-DMD modes have a slightly higher energy (−10.9dB and −12.9 dB) compared to mode two (−11.1 dB and −14.0 dB).The other cases do not show this behavior.

Opt-DMD Coherent Structures
In addition to characterizing oscillation frequencies and growth or decay rates, DMD also provides insight into wake structure.Reconstructions for the first two modes are shown in figure 8.Each modal reconstruction corresponds to the phase of oscillation shown in figure 5.For all St values, mode 1 exhibits top-bottom symmetry in vortex structure across the centerline (y/c = 0).In contrast, mode 2 has an antisymmetric structure with vortex lobes of smaller size compared to mode 1.The higher order modes 3 and 4 are also symmetric and asymmetric respectively, as shown in figure 9. Qualitatively, these modes better capture skews in the wake structure across the centerline, particularly the for the case of St = 0.16 (figures 9a-b).It was shown by Moored et al. (2012) that the overlap in vorticity modes at certain locations in space contributes to producing the full vortex structure.Similar observations can also be made in the present study, where the modal hierarchy is trying to represent a convective process in terms of purely oscillatory DMD modes that are harmonics of the flapping frequency.This is reflected in the DMD mode structure: since mode 2 exhibits twice the oscillation frequency as mode 1, vortex lobes from mode 2 exhibit roughly twice the spatial wavenumber as mode 1 (see figure 8).A coherent vortex is formed at a particular region in space when both modes are either in-phase or in anti-phase and vortex lobes with the same sign overlaying each other.This is also consistent with the fact that the foil sheds a vortex each half cycle.Whether the vortex formed is positive or negative is determined by the sign of the spatial eigenfunctions.These features (and the modes themselves) are phase-locked in a reference frame moving with the local convective velocity.
Notable differences between the modes from the near-momentumless wake and that from the 2P and 2S wakes are observed.For instance, mode 1 for St = 0.16 has approximately four symmetric lobes of vorticity along the lateral axis (figure 8a).These elongated vortex packets expand laterally as they advect downstream and similar trends can be seen for mode 2 (figure 8b).This is consistent with the lateral variation observed in the mean velocity profile and phase-averaged vorticity fields.In contrast, mode 1 for the cases St = 0.23, 0.29, and 0.35 primarily consists of two vortex lobes with decreasing streamwise extent.Thus, the size of the lobes reflects the length scale of the shed vortices, which is partially determined by the convective length scale U ∞ / f .The differences that reflect the transition from the 2P wake to the 2S reverse von Kármán street are more subtle.In the case of St = 0.23 which exhibits a 2P wake morphology, the vorticity lobes travel consistently downstream (see figures 8c, d).As the Strouhal number increases from St = 0.23 to St = 0.35, the lobe pairs move inward towards the centerline over the region 0.22 ≤ x/c ≤ 1.0.In the region where x/c > 1.0, the lobes consistently travel in the streamwise direction again.This dynamic characteristic can be observed better from modes 3 and 4 in figure 9(c − h), where the lobes structures are considerably smaller.Another key distinction across the 2P to 2S wake transition is that the symmetric vorticity lobes of modes 1 and 3 start to coalesce into a single lobe towards the downstream end of the field of view, x/c ≥ 2.0.The antisymmetric pairs for modes 2 and 4 remain apart.This is most noticeable for the case of St = 0.35 (figures 8g and 9g).
Although the high-frequency wake structures exhibit reverse von Kármán streets -as expected for thrust-producing systems -their symmetric and antisymmetric topologies are similar to that from the classical Kármán vortex street shed from twodimensional objects (Tu et al., 2014;Taira et al., 2020;Araya et al., 2017).In these drag-producing wakes, the dominant DMD or POD modes obtained for these flows can also correspond to a series of frequency harmonics.One distinction from the classical drag-producing bluff body wakes is that their symmetric mode structures may consist of a single lobe of vorticity across the y-axis -as shown in experiments from (Tu et al., 2014) -rather than two or more from our study.

Influence of Modes on Wake Dynamics
To further understand how the symmetric and antisymmetric modes relate to the full flow field, we superimpose the streamwise component associated with DMD modes 1 ( ũ1 ) and 2 ( ũ2 ) separately with the time-averaged streamwise velocity ( ū).Figures 10 and 11 show the full opt-DMD reconstructions (i.e., comprising modes 1-4) as well as these single-mode reconstructions of the streamwise velocity component.In the original streamwise fields, the vortex development is characterized by the regions where u is less than the that of the freestream velocity (u/U ∞ < 1), resulting in shear layer roll up.In the transition to the 2S wake, a wavy jet is observed in the original flow fields (11a, b).Mode 1 with the mean flow (u = ū + ũ1 ) reproduces most of these wake dynamics, including a reasonable portion of the shear layer roll up from which The wake dynamics associated with mode 1 and the associated symmetric vorticity perturbations closely resemble the spatial instabilities observed by Moored et al. (2012) in the 2S wake structure created by a flexible fin.The main distinction lies in the number of vorticity lobes, with the former exhibiting two lobes instead of three.This observation suggests that the opt-DMD modes yield coherent structures that are related to spatial instabilities.It is also likely that mode 1 induces a majority of the net thrust since this mode is associated with the lower velocity structures that form the vortex structures outside the central jet.In contrast, mode 2 is associated with the shedding frequency of each vortex or vortex pair for the 2S and 2P wakes (i.e., twice the oscillation frequency).Mode 2 accounts for much of the remaining shear in the center jet region; this is particularly evident in the higher Strouhal number wakes (figure 11e − f ).

Coherent Reynolds Stress Contributions from DMD Modes
The subtle differences in DMD mode structure across Strouhal numbers suggest that the associated Reynolds shear stresses may provide additional insight into thrust and drag effects.Following Reynolds and Hussain (1972) and assuming that the Reynolds stress contributions from the turbulent fluctuations u ′ are negligible compared to the contributions from the phase-averaged periodic components ũ, the time-averaged momentum equation in the streamwise direction can be expressed as follows: The equation above also assumes no external pressure gradient, a purely 2D flow, and that cross-stream gradients of the viscous and Reynolds stress terms dominate over the streamwise gradients.Invoking both a strong parallel flow assumption, whereby v ≈ 0 and ∂ ū/∂x ≈ 0 locally, equation 18 can be further simplified to: which can be used to estimate the induced mean flow.In the above equation, ūn can be thought of as the mean velocity induced by DMD mode n, i.e., the Reynolds shear stresses generated by the periodic velocity components ũn and ṽn associated with mode n.Since the DMD modes are harmonics of the foil oscillation frequency, there are no mean Reynolds shear stress contributions from interactions across modes.In other words, we expect ũm ṽn = 0 for m ̸ = n.Equation ( 19) indicates that positive values of d ũn ṽn dy are associated with minima in ūn and vice versa.Positive ūn contributions increase momentum in the wake and produce thrust.The opposite is true for negative ūn contributions, which result in induced drag.
Predictions made using the highly simplified momentum balance in equation 19 are complementary to the control volume analyses of momentum entrainment and expulsion presented by Moored et al. (2014), which showed that the surrounding fluid is entrained into the near wake region close the foil at wake resonance.The increase in mass flow rate of the jet wake produces thrust, and the entrainment region should therefore have a Reynolds stress distribution in which the periodic components ũ ṽ induce a jet-like mean profile.We solve equation 19 numerically for ūn using the Reynolds shear stress profiles for DMD modes 1 and 2 that are closest to the foil (x/c = 0.22), enforcing the boundary conditions ūn (−∞) = ūn (∞) = 0.
Figures 12 and 13 respectively show the coherent Reynolds stress fields ũn ṽn and induced mean profiles ūn predicted using equation ( 19) for opt-DMD modes 1 and 2. The Reynolds stress fields show patterns that are antisymmetric across the centerline over the region of 0.22 < x/c < 1.0.This antisymmetry begins to break down for the case of St = 0.35 (figure 12g and 13g).As expected, mode 1 (i.e., the primary harmonic) generates larger Reynolds shear stresses and induced mean flow contributions compared to mode 2. However, this discrepancy is more pronounced for the higher Strouhal number cases (St = 0.29, 0.35) exhibiting 2S-type wakes.The maximum magnitudes of ũ1 ṽ1 for mode 1 increase with Strouhal number.
The shape of the induced mean flow profiles for mode 1 are similar for the St = 0.16 and 0.23 cases (figure 12b, d), where two distinct positive maxima can be seen (with one noticeable negative component for the case St = 0.23).Induced mean profiles for both 2S wakes (St = 0.29, 0.35) in 12( f , h) also show similar characteristics.However, ū1 for St = 0.29 shows two separate maxima while the profile for St = 0.35 shows a single maximum near the centerline.This transition from two distinct maxima in the mean profile to a single peak is consistent with the 2P-2S transition noted earlier as the Strouhal number increases from St = 0.23 to St = 0.35.Interestingly, ū1 profiles for St = 0.29 and St = 0.35 match the measured mean profile shapes shown in figure 5(c, d).Specifically, the St = 0.29 profile is approximately trapezoidal in shape while the St = 0.35 profile resembles a typical jet.Induced ū1 profiles for these high Strouhal number cases are also much larger in magnitude compared to the lower Strouhal number cases (St = 0.16, 0.23), which is indicative of higher entrainment and thrust production.
Generally, Reynolds shear stress contributions and induced velocities from mode 2 in figure 13 are lower than those from mode 1. Particularly, ū2 ≪ ū1 for the 2S wakes in St = 0.29 and St = 0.35.Additionally, While the mode 1 contributions vary significantly in magnitude across Strouhal number, all of the mode 2 profiles exhibit similar maximum values of ũ2 ṽ2 and ū2 .Reynolds stress fields associated with mode 2 for St = 0.16 and St = 0.23 (figure 13b, d) remain reasonably coherent over the field of view.The resulting mean flow profiles also show consistent negative values ( ū2 < 0), which is indicative of drag generation.This aligns with Floryan et al. (2017), who found that scaling laws differed from experiments in efficiency due to viscous drag effects at lower Strouhal numbers.the coherent stress fields are likely characteristics of the transition to bluff body shedding.Note that these drag- inducing effects are characterized by modes with a higher frequency than that from typical bluff body shedding (c.f., similar observations from Strom et al., 2022).
While the ū2 profiles are primarily negative for the lower Strouhal number cases, distinct positive regions are observed as the flow transitions from the 2P wake to a 2S wake for St = 0.29 and 0.35 as seen in figure 13( f , h).For St = 0.29, the ū2 profile is positive and approximately trapezoidal in shape, matching the mean flow in figure 5 f .For St = 0.35, the induced ū2 profile has a mix of positive and negative regions, with the total area under the profile yielding a net negative value.The transition from purely negative ū2 contributions for St = 0.16 and 0.23 to primarily positive ū2 contributions for St = 0.29 is reminiscent of the earlier convective instability workas noted by Triantafyllou et al. (1993), in a reverse von Kármán vortex street, there is little to no competition between the drag wake and the thrust wake.This is consistent with the present observations for St = 0.29.

Deterioration of Propulsive Efficiency at High St
The marginal deterioration in propulsive efficiency along with the negative ū2 contributions for the St = 0.35 case can potentially be attributed to effects from leadingedge vortices.While we do not have access to PIV data around the foil, force and torque measurements, along with the angle of attack trends support this interpretation.In particular, we expect these effects to be evident in the normal or lift forces on the foil, which also influence the power input and propulsive efficiency (see equations 4 and 5).Dimensionless phase-averaged lift forces F y and power requirements ℘ measured for Strouhal numbers St = 0.29 − 0.59 are shown in figure 14.As expected, the lift forces exhibit the same oscillation frequency as the actuation, which also matches the first mode frequencies.The relative power exhibits oscillations at twice the actuation frequency, which is also the second mode frequency.
Noticeable deviations in relative lift and power are observed between the nearoptimal Strouhal numbers (St = 0.29 and 0.35) and the highest Strouhal number case (St = 0.59).Particularly, the small enhancements in the magnitude of the lift force observed between t/T ≈ 0.4 − 0.5 and t/T ≈ 0.9 − 1.0 lead to lower relative power requirements, as shown in figures 14(b, c).For the highest propulsive efficiency case of St = 0.29, additional decreases in relative power at t/T ≈ 0.25 and 0.75 lower the overall power across the full oscillation cycle (figure 14c).Note that the relative lift force enhancements for this Strouhal number case appear soon after maxima and minima in the effective angle of attack (t/T = 0.25, 0.75).
These relative lift force enhancements may be linked to the structure of the Reynolds shear stress contribution and positive induced velocity from DMD mode 2 for St = 0.29 (figures 13e, f ), suggesting a delay in separation of the leading edge vortices.The separation effects may start to occur for the St = 0.35 case, which show a negative induced velocity from mode 2 (figure 12h) and higher relative power requirements in comparison to those from St = 0.29, resulting in lower propulsive efficiency.The deterioration in performance is even more visible for the St = 0.59 case, which shows a substantial decrease in the magnitude of F y /F y0 and corresponding increases in ℘/℘ 0 .
The tapering of efficiency past St = 0.29 can also be explained based on high angles of attack.The stall angle for a static NACA 0012 foil is approximately 16 • .However, oscillating loads can increase the stall angle substantially (Maresca et al., 1979).For the present experiments, the maximum angles of attack α 0 for the higher Strouhal number cases St = 0.29, 0.35, and 0.59 are 28.2 • , 33.1 • , and 47.5 • , respectively.Thus, it is likely that the induced drag effects from the secondary mode (figure 13h, f ) become more relevant as the effective angle of attack increases, through which separation effects may occur -as also evident in the F y /F y0 traces.In turn, the efficiency tapers downward past the Strouhal number St = 0.29.However, it would be misleading to say that all cases with high angles of attack lead to reduced efficiency, as some combinations of kinematic parameters can delay these dynamic stall effects (Anderson et al., 1998;Maresca et al., 1979;Ellington et al., 1996).Instead, a high angle of attack where the flow stays attached to the foil can potentially maximize efficiency, an effect that has also been observed for flexible propulsors (Quinn et al., 2015).

Discussion and Conclusions
A triple decomposition method in which the periodic component of the wake is composed of opt-DMD modes was used to provide further physical insight into the propulsive performance of oscillating foils.The experimental data in which the method was used are broadly consistent with prior literature.We observe a near-momentumless wake structure for St = 0.16 and peak propulsive efficiency for St = 0.29.Over this range of Strouhal numbers, PIV measurements show a transition from a momentumless wake to a 2S wake morphology associated with the classical reverse von Kármán vortex street.The opt-DMD method is employed here instead of the classical or (exact) DMD approach as it captures more energy in the periodic component of the flow field.For this particular foil, the opt-DMD modes appear as harmonics of the foil oscillating frequency, with alternating symmetric and antisymmetric morphology across the wake centerline.
The most interesting finding from our study is the relationship between the coherent Reynolds stress contributions from DMD modes and their impact on propulsive performance.Although equation 19 is only an approximation to the mean flow equation for the wake, it helps delineate modal contributions to drag and thrust.In our experiments, the mean velocities induced by the primary opt-DMD mode (mode 1) generally show thrust-producing characteristics across all St values, though there is a transition from a two-hump mean profile for the 2P wakes to a jet-like mean profile for the 2S wakes.However, several different trends emerge from the mean velocities induced by mode 2, representing features oscillating at twice the fundamental frequency of the foil: (i) for Strouhal numbers below the optimum (St < 0.29) the induced mean velocities ū2 are negative and suggest a transition to bluff body shedding; (ii) the secondary mode for the optimum Strouhal number case (St = 0.29) has a prominent positive induced mean velocity, indicating that the thrust wake is dominant; (iii) past the optimum Strouhal number, mode 2 transitions back to a negative induced mean velocity, which may be indicative of flow separation as a result of high angles of attack (St > 0.29).
Given the close connection between DMD and stability analyses, it may be of interest to compare the opt-DMD modes with features identified in prior stability analyses (Moored et al., 2012;Arbie et al., 2016).The development of the opt-DMD modes close to the foil suggest that their may be some spatial growth of these perturbations in the near-field.Thus, there may be links between these opt-DMD modes and the spatially-growing modes identified in the 'wake resonance' studies (Triantafyllou et al., 1993;Lewin and Haj-Hariri, 2003;Moored et al., 2014).Nonetheless, it should be emphasized that the temporal DMD method was used in this study, which extracts modes at a particular oscillation frequency.A spatial DMD approach could provide a more direct approach to understanding these instability mechanisms.In this case, a spatial resolution that is higher than the the PIV data collected in this study would be necessary to allow for an appropriate spatial DMD analysis.
This study focused on a narrow parameter range: we studied a single, rigid foil exhibiting periodic single-frequency oscillations.However, modal analysis methods similar to those employed here could provide substantial insights into wakes involving more complex kinematics, dynamic fluid-structure interactions, multiple-foil interactions, or massively separated flows (see e.g., Raspa et al., 2013;Andersen et al., 2017).

Figure 1 :
Figure 1: Geometry and parameters for a pitching and heaving hydrofoil.

Figure 4 :
Figure 4: Comparison of time-averaged thrust coefficients and propulsive efficiencies as a function of Strouhal number.

Figure 7 :
Figure 7: Power spectral densities for the original flow field (black) and the reconstructed flow fields from the opt-DMD (orange) and exact-DMD (blue) method for (a) St = 0.16, (b) St = 0.23, (c) St = 0.29, and (d) St = 0.35.The dotted lines represent exponential decay rates obtained from fits to the first four mode amplitudes from opt-DMD and show the following behavior: e −5.1St (a), e −4.0St (b), e −1.5St (c), and e −1.0St (d).

Figure 14 :
Figure 14: Dimensionless heave cycle h/c (dotted line) and effective angle of attack α (solid lines) for one cycle (a), together with phase-averaged lift forces F y (b) and power requirements ℘ (c) normalized by maximum values (F y 0 and ℘ 0 ) for St = 0.29 (red), St = 0.35 (blue), and St = 0.59 (black).

Table 1 :
Parameters used in the present and previous oscillating-foil studies.Chordbased Reynolds number Re c = ρU ∞ c/µ, Strouhal Number, St, phase between heaving and pitching ϕ p , pitch amplitude θ 0 , and dimensionless heave ratio h * = h 0 /c.