Near-wake dynamics of a vertical-axis turbine

Abstract

Cross-flow, or vertical-axis, turbines are a promising technology for capturing kinetic energy in wind or flowing water and their inherently unsteady fluid mechanics present unique opportunities for control optimization of individual rotors or arrays. To explore the potential for beneficial interactions between turbines in an array, as well as to characterize important cycle-to-cycle variations, coherent structures in the wake of a single two-bladed cross-flow turbine are examined using planar stereo particle image velocimetry in a water channel experiment. There are three main objectives in the present work. First, the mean wake structure of this high chord-to-radius ratio rotor is described, compared with previous studies, and a simple explanation for observed wake deflection is presented. Second, the unsteady flow is then analysed via the triple decomposition, with the periodic component extracted using a combination of traditional techniques and a novel implementation of the optimized dynamic mode decomposition. The latter method is shown to outperform conditional averaging and Fourier methods, as well as uncover frequencies suggesting a transition to bluff-body shedding in the far wake. Third, vorticity and finite-time Lyapunov exponents are then employed to further analyse the oscillatory wake component. Vortex streets on both sides of the wake are identified, and their formation mechanisms and effects on the mean flow are discussed. Strong axial (vertical) flow is observed in vortical structures shed on the retreating side of the rotor where the blades travel downstream. Time-resolved tracking of these vortices is performed, which demonstrates that vortex trajectories have significant rotation-to-rotation variation within one diameter downstream. This variability suggests it would be challenging to harness or avoid such structures at greater downstream distances.

1. Introduction

Cross-flow, or vertical-axis, turbines are experiencing a resurgence in research interest for the conversion of wind and water currents to electricity. One motivation is the mounting evidence that arrays of closely spaced cross-flow turbines can extract more energy per unit land area than industry-standard axial-flow turbines. This property has benefits where the array mounting area is limited, such as roof-top installations, and where the region of high flow speed is concentrated, such as mountain passes or tidal channels. Dabiri (2014) and Dabiri et al. (2015) reported a power density of 10–20 W m$^{-2}$, compared with the 1–3 W m$^{-2}$ output of conventional axial-flow wind turbine farms (MacKay 2008; Adams & Keith 2013). Notably, in field experiments, Brownstein, Kinzel & Dabiri (2016) find average rotor performance in an array is 20 % higher than the performance of a single isolated turbine. Similarly, Scherl et al. (2020) have demonstrated a 30 % increase in average rotor output in an array of two cross-flow turbines compared with isolated turbine performance. The remarkable performance of densely packed cross-flow turbines stems from several fluid mechanical phenomena. First, the orientation of the rotation axis results in an acceleration of the bypass flow, especially on the side of the rotor where the blades are retreating (travelling downstream). Neighbouring rotors placed in this flow benefit from the increased incident mean velocity. Second, the tip vortices shed from the blades have an axis of rotation that lie in a plane parallel to the ground. These vortices induce vertical mixing, increasing the transfer of momentum from the high-speed flow above the array to the rotor level, increasing the streamwise wake recovery rate (Bachant & Wosnik 2015). Finally, we speculate that performance may be enhanced through the interaction between periodic coherent structures shed by an upstream turbine and the blades of a downstream turbine. This hypothesis is inspired by schooling fish who have been shown to benefit from well-timed interactions with vortices shed from upstream individuals (Whittlesey, Liska & Dabiri 2010; Maertens, Gao & Triantafyllou 2017). These potential performance increase mechanisms motivate the study of the mean and periodic components of a cross-flow turbine wake, with a special focus on coherent structures that may interact with nearby turbines in an array and be exploitable through control.

Measurement and analysis of cross-flow turbine wakes have been conducted for decades, starting with Muraca & Guillotte (1976). Point measurements using Pitot tubes (Muraca & Guillotte 1976; Battisti et al. 2011), hot-wires (Bergeles, Michos & Athanassiadis 1991; Battisti et al. 2011; Peng, Lam & Lee 2016; Persico et al. 2016), laser (Buchner et al. 2018), acoustic (Kinzel, Mulligan & Dabiri 2012; Bachant & Wosnik 2015; Kinzel, Araya & Dabiri 2015) and Doppler velocimetry have been used to describe the mean wake structure, spectra and time-average turbulence statistics. Two-component (Araya & Dabiri 2015; Eboibi, Danao & Howell 2016; Posa et al. 2016; Araya, Colonius & Dabiri 2017) and three-component (Tescione et al. 2014; Rolin & Porté-Agel 2015; Hohman, Martinelli & Smits 2018) planar particle image velocimetry (PIV) and magnetic resonance velocimetry (Ryan et al. 2016) measurements, as well as simulations (Scheurich, Fletcher & Brown 2011; Scheurich & Brown 2013; Nini et al. 2014; Shamsoddin & Porté-Agel 2014; Boudreau & Dumas 2017), have been used to investigate the wake spatial variability, including wake geometry, recovery rate and the role of turbulence. Despite the widely varying rotor configurations and operating conditions across these studies, a set of features common to cross-flow turbine wakes have emerged. First, wake measurements have often been made at mid-plane of the rotor, perpendicular to the rotation axis. In this plane, all studies report some asymmetry or angular deflection of the wake in the direction of turbine rotation, with a more intense shear layer on side where the blades are advancing (travelling upstream). Flow structures shed at the blade passing frequency (rotational frequency $x$ number of blades) have been identified in nearly all studies on the retreating side of the wake (e.g. Brochier et al. 1986; Battisti et al. 2011; Bachant & Wosnik 2015; Ryan et al. 2016) and on both the retreating and advancing sides (e.g. Posa et al. 2016; Araya et al. 2017; Boudreau & Dumas 2017; Hohman et al. 2018). Araya et al. (2017) determined that these shear flow oscillations transition to those corresponding to a bluff body in the far wake. Second, areas of high turbulence intensity have been identified in a streak on the advancing side of the wake deficit (Bachant & Wosnik 2015) and on both sides (Rolin & Porté-Agel 2015; Hohman et al. 2018). Third, studies that examined the three-dimensional wake structures have identified the primary mechanism for wake recovery as axial (vertical) flow induced by vortices shed from the blade tips (Kinzel et al. 2012, 2015; Boudreau & Dumas 2017) or the induced cross-stream (horizontal) flow (Bachant & Wosnik 2015). In contrast, the wake recovery in axial-flow turbines is driven primarily by turbulent mixing upon the breakdown of the helical tip vortices (Lignarolo et al. 2015; Boudreau & Dumas 2017). Consequently, wake recovery rates have been documented to be significantly faster than those of axial-flow turbines (Dabiri 2011; Boudreau & Dumas 2017). In addition to wake measurements, multiple studies have performed measurements within the rotor, demonstrating the importance of dynamic stall and subsequent blade–vortex interactions in normal cross-flow turbine operation (Brochier et al. 1986; Fujisawa & Shibuya 2001; Ferreira et al. 2009; Edwards, Danao & Howell 2015; Eboibi et al. 2016; Dave et al. 2021).

The contributions of the present study are threefold. First, time-average, three- component, planar PIV measurements are presented for a relatively high chord-to-radius ratio turbine and the wake structure is related to the turbine rotor hydrodynamics. Second, we demonstrate that an algorithm incorporating the dynamic mode decomposition (DMD) can identify energetically important modes that cannot be discovered by other methods. Third, by analysing the form and trajectory of coherent structures shed into the near wake, we identify the region over which their propagation is deterministic, which is of relevance to array control.

2. Cross-flow turbine operation and experimental methods

2.1. Cross-flow turbine background

Despite typically having only a single degree of freedom, rotation about a central axis, the fluid dynamics of cross-flow turbines is inherently unsteady. This is because, even with a steady inflow, the local flow conditions experienced by the blade vary cyclically over the course of a single rotation. Neglecting flow variations induced by the turbine, the local flow velocity magnitude and angle of attack vary according to

(2.1)\begin{equation} U_n(\theta)^* = \frac{| \boldsymbol{U_n}(\theta) |}{U_\infty} = \sqrt{\lambda^2+2\lambda \cos(\theta)+1}, \end{equation}

and

(2.2)\begin{equation} \alpha_n(\theta) ={-}\text{Tan}^{{-}1}[\sin(\theta), \lambda + \cos(\theta)]+\alpha_p, \end{equation}

respectively, where $\text {Tan}^{-1}$ is the four quadrant arctangent and $\alpha _p$ is the preset pitch (blade mounting) angle, $\theta$ is the blade azimuthal position and $\lambda$ is the tip-speed ratio, or non-dimensional rotation rate; $\lambda$ is given by

(2.3)\begin{equation} \lambda = \frac{\omega R}{U_\infty},\end{equation}

where $\omega$ is the angular velocity of the turbine, $R$ is the radius and $U_\infty$ is the free-stream velocity. Although it is possible to approximate the effect of induction as in Ayati et al. (2019), the authors are not aware of an approximation that includes both streamwise and cross-stream induction. Cross-stream force analysis presented in § 3 indicates significant induction that should not be ignored in an induction approximation. As shown in figure 1, we define $\theta = 0$ where the quarter-chord of the blade is travelling directly upstream. Vector diagrams of these quantities, as well as examples of how they vary over the course of one rotation are given in figure 1. Especially at low tip-speed ratios, the local angle of attack can far exceed the static stall angle of the foil. These kinematics, equivalent to a rapid pitching manoeuvre, can lead to dynamic stall and the corresponding roll-up of a leading edge vortex (Eldredge & Jones 2019). Depending on the timing, strength and trajectory, this vortex may contribute to or detract from power output (Ferreira et al. 2009; Strom, Brunton & Polagye 2015). Dynamic stall and the resulting coherent structures provide an opportunity for optimization of blade–fluid structure interactions, either in the case of a single turbine, as in Strom, Brunton & Polagye (2017), or for multi-rotor interactions. Array optimization that seeks to maximize power based on the mean flow has been successfully demonstrated (Brownstein et al. 2016). However, arrays of cross-flow turbines may also be able to take advantage of periodic coherent structures in the wakes of nearby turbines (Scherl et al. 2020). In addition to their potential influence on array interactions, the role of coherent wake structures in deficit recovery rate motivates the close examination of their lifetime and trajectory.

Figure 1. (a) Diagram of geometric and kinematic quantities. (b) Free-stream, rotational and resulting total velocity vector and local angle of attack. (c) Variation in local angle of attack (top) and flow velocity (bottom) as a function of azimuthal blade position for three values of tip-speed ratio ($\lambda$).

The efficiency with which a cross-flow turbine converts flow kinetic energy to rotational mechanical energy is given by

(2.4)\begin{equation} C_P = \frac{\omega \tau}{\frac{1}{2} \rho U_\infty^3 A}, \end{equation}

where $\omega$ is the turbine rotation rate, $\tau$ is the mechanical torque produced by the rotor, $\rho$ is the operating fluid density and $A$ is the rotor swept area. Performance is often characterized as a function of the tip-speed ratio in (2.3) and can be presented as a time-average or phase-averaged quantity (Polagye et al. 2019). For wake measurements to be most meaningful, it is useful to operate the turbine at a realistic operating condition, such as the tip-speed ratio that yields the maximum $C_P$, because wake characteristics differ significantly between power producing and non-power producing operating states (Araya & Dabiri 2015).

At large turbine scales, such as those used for commercial power production, rotor geometries with few, relatively small chord-length blades exhibit high maximum efficiency. One example is the Sandia 34 m test bed turbine, with a peak $C_P$ of approximately 0.41 (Ashwill 1992). However, because cross-flow turbine performance can improve rapidly with increasing Reynolds number (Bachant & Wosnik 2016; Miller et al. 2018), at smaller scales it can be useful to increase the chord length of the turbine to maximize the blade Reynolds number. Larger chord-length foils may also be more structurally robust, which is important because blade fatigue is often the cause of cross-flow turbine structural failure (Möllerström et al. 2019). Finally, large chord-to-radius turbines have a lower tip-speed ratio at peak performance, reducing losses from support structures and radiated noise. These factors motivate the study of cross-flow turbines with relatively high chord-to-radius ratios. The low tip-speed ratio at peak efficiency of this geometry results in large local angle-of-attack variations (figure 1) and in separation and stall when operating at maximum $C_P$ (Snortland, Polagye & Williams 2019).

2.2. Flume and turbine

Experiments were performed in the Alice C. Tyler flume at the University of Washington. The flume has a test section measuring 0.76 m wide and 4.9 m long. The dynamic water depth was 0.47 m, the free-stream velocity ($U_{\infty }$) was maintained at 0.7 m s$^{-1}$ and the turbulence intensity (${rms}(\boldsymbol {u}')/U_\infty$) was 1.5%. The water temperature was held at a constant 16.3$\pm$0.4 $^\circ$C.

The cross-flow turbine model had a height of $H = 0.234$ m and a diameter of $D = 0.172$ m. The diameter-based Reynolds number,

(2.5)\begin{equation} Re_D = \frac{D U_\infty}{\nu}, \end{equation}

was $1.1\times 10^5$, where $\nu$ is the water kinematic viscosity. The turbine was vertically centred in the flume with a depth-based Froude number

(2.6)\begin{equation} Fr = \frac{U_\infty}{\sqrt{g d}}, \end{equation}

was 0.33 where $g$ is the acceleration due to gravity and $d$ is the dynamic water depth. The blockage ratio, or the ratio of the turbine cross-sectional area to the test section cross-sectional area, was 11 %. The turbine consisted of two, straight, NACA0018 profile blades with chord length of $c=0.061$ m for a chord-to-radius ratio of 0.71 and solidity,

(2.7)\begin{equation} \sigma = \frac{N c}{{\rm \pi} D}, \end{equation}

of 0.225, where $N$ is the blade count. The blades were mounted to a 0.012 m diameter central shaft via circular endplates at a pitch angle of 6$^\circ$ (leading edge rotated outwards about the quarter-chord). The turbine was operated under constant angular velocity control at its peak performance point of $C_P = 0.26$ at a tip-speed ratio of $\lambda = 1.2$. The turbine performance curve is given in figure 2, and details on the methods used to determine turbine performance can be found in Strom et al. (2017).

Figure 2. Performance curve (mechanical efficiency vs tip-speed ratio) for the experimental turbine. Red cross indicates operating point during wake data collection. Performance curve data were collected using the experimental set-up detailed in Strom et al. (2017), but during the PIV experiments with a cantilevered turbine (figure 3), the upper load cell was removed to increase the stiffness of the experimental set-up.

The turbine was cantilevered from the face of a direct-mount servomotor (Yaskawa SGMCS) with an integrated 1048576 edges-per-revolution encoder providing blade position feedback, which was recorded via a counter on a National Instruments PCIe data acquisition card to a computer at a rate of 1 kHz. This acquisition was synchronized with the PIV measurement system. The servomotor regulated the rotational speed of the turbine to a constant value and power generated was actualized as reverse current in the servomotor and dissipated in a dump resistor. The turbine rotor and servomotor were mounted to a robotic gantry system, providing accurate translation of the rotor in the streamwise direction.

In a separate set of experiments described in Strom et al. (2017), two six-axis load cells are used to measure the reaction forces between the turbine system and both the flume floor and an upper mounting beam. As a result, streamwise ($F_X$) and cross-stream ($F_Y$) forces are measured, and the corresponding coefficients are calculated as

(2.8a,b)\begin{equation} C_{FX} = \frac{F_X}{\frac{1}{2}\rho U_\infty^2 A} \quad \text{and} \quad C_{FY} = \frac{F_Y}{\frac{1}{2}\rho U_\infty^2 A}, \end{equation}

respectively.

2.3. PIV measurement

Measurements of the turbine wake were acquired using time-resolved stereo planar PIV. Data were collected in a free-running manner at 100 Hz, corresponding to 5.35$^\circ$ of blade rotation between measurements, and was not locked to specific blade positions. Measurements were taken at the mid-span of the turbine, in the plane normal to the axis of rotation. Illumination was provided by a Continuum TerraPIV Nd:YLF laser, and images were captured by two Phantom V641 cameras, with resolutions of $2560 \times 1600$ pixels. Cavitation bubbles from the flume recirculation pump were used as passive tracers and measured approximately 1.5 pixels in diameter. Measurement resolution was increased by capturing the wake using six overlapping fields of view, as illustrated in figure 3(a). We note that data are collected in these six regions in separate experiments, and therefore not synchronized. Consequently, we align the data in post-processing using the algorithm discussed in § 4 and presented in more detail in Nair et al. (2020). The combined measurement area, shown in figure 3(b), started 0.57$D$ downstream from the turbine axis, and extended 3.68$D$ downstream, and 3$D$ in the cross-stream direction.

Figure 3. Turbine and PIV measurement set-up diagram (a) and PIV measurement locations in the mid-plane along the $\boldsymbol {z}$ direction (b).

Spatial calibration was performed with custom stereo calibration target spanning the entire width of the flume section in conjunction with a robotic camera gantry used to repeatably move the cameras in the cross-stream direction. Post-processing was performed with custom image manipulation software and TSI Insight for the cross-correlation. Ghost velocities due to small laser-sheet/calibration target misalignment was corrected through image warping in post-processing. Velocity fields were calculated using iterative multi-grid processing, with initial square interrogation window side size of 64 pixels and a final size of 16 pixels. With 50 % window overlap, the resulting velocity vector spacing was 0.0068$D$.

In addition to the mid-plane fields that are presented in this work, velocity fields were also collected above and below the mid-plane, as shown for vertical velocity in figure 4. A visual comparison of the three planes shows a strong vertical flow and asymmetry in the wake. Although this is an interesting observation, the present work is restricted to the mid-plane for two key reasons. First, this work primarily serves as a methodological exploration of the mean, phase-averaged and time-resolved wake structure. Expanding this to encompass multiple vertical planes would be unwieldy. Second, and more importantly, this asymmetry appears to originate within the confines of the rotor. As such, the field of view captured in this set of experiments cannot definitively identify the origin of this asymmetry – it can only confirm the presence of strong vertical flow and asymmetry in the wake. Collecting in-rotor flow fields may answer this question, and this is the subject of ongoing work.

Figure 4. Vertical variation in vertical velocity (arbitrary phase).

3. Results and analysis: mean flow

The mean wake deficit contours and normalized velocities are shown in figure 5. As in prior work, we observe an asymmetric wake deficit with an intense shear layer on the advancing side of the wake (see figure 3(b) for advancing vs retreating nomenclature). Wake deficit recovery occurs faster on the retreating side, as previously observed by Tescione et al. (2014). The mean wake deficit is never negative, meaning there is no recirculation region. Araya et al. (2017) showed a decrease in wake deficit with reduction in the number of blades; however, even their two-bladed turbine showed some negative streamwise velocity. A survey of wake measurements in prior work indicates that neither rotor efficiency, solidity nor the expression of dynamic solidity of Araya et al. (2017) are good predictors of whether or not a negative wake deficit occurs. It is possible that some combination of these factors, in addition to the test section blockage ratio and rotor thrust, would be necessary to predict the magnitude of the wake deficit.

Figure 5. (a) Mean wake deficit profiles. Streamwise velocity profiles along cross-stream stations (dashed lines). The distance from one station to the next is a change in velocity equivalent to the mean free-stream velocity, $U_\infty$. (b) The mean streamwise, (c) cross-stream and (d) vertical (axial) velocities, normalized by the free-stream velocity.

Despite the differences in turbine geometry, the streamwise wake velocity is similar to those described by Peng et al. (2016) (five blades, $c/R = 0.3$) and Hohman et al. (2018) (three blades, $c/R = 0.2$), suggesting that rotor geometry has limited effect on the mean wake structure. The blockage ratio of 11 % also increases the shear between the wake and bypass flow compared with an unconfined case, but, as noted in Ross & Polagye (2020), the time-average wake structure is qualitatively invariant with this magnitude of blockage.

There are several conflicting theories in the literature about the root cause of the mean wake profile asymmetry. Specifically:

• • Araya et al. (2017): ‘In all cases, there is a notable asymmetry of the [cross-flow turbine] wake. This is attributed to the stronger shear layer that forms on the side of the turbine where the blades are advancing upstream’.

• • Hohman et al. (2018): ‘This behavior is as expected, as the majority of the power is generated on the advancing side of the turbine, and therefore a larger momentum deficit will be seen on this side’.

• • Bachant & Wosnik (2015): compares the effect with that of a rotating cylinder, stating ‘Compared with the rotating cylinder wake measurements of Lam & Peng (2016) we see a similar asymmetry in the mean streamwise velocity. The wake is less asymmetrical with respect to the wake centreline for the turbine compared to the rotating cylinder for the same non-dimensional rotation rate, although some of these differences may be due to the cylinder experiments’ lower Reynolds numbers’.

• • Peng et al. (2016): provide multiple explanations: ‘There are two major factors that may contribute to this wake asymmetry. One factor is that more turbulent structures are produced at the windward than at the leeward. When the blade advances under adverse pressure gradients at the windward, stronger vortex shedding and much severer flow separations take place. The other factor is that the wake flows are transported toward the windward. First, when the blade moves upwind at the windward, it causes stronger blockage effect compared to that at the leeward. Therefore, at the windward, the blade wake is characterized by a lower pressure, which induces the cross-wind flows. Second, when the blade operates at the downstream half-revolution, the strong angular momentum drags and propels the wake flows toward the windward’.

We propose three mechanisms for the wake asymmetry. The first two revolve around the balance between blade forcing and flow acceleration governed by Newton's second law. The third involves shed vorticity and is discussed in § 5.2. The first mechanism producing asymmetry is a difference in streamwise forcing between the advancing and retreating sides of the rotor. On average, the streamwise forcing on a blade is greater on the advancing side, see figure 6(a) (yellow vectors) and (b). This is due to a combination of the difference in direction and magnitude of the relative translation of the blade compared with the free-stream flow (compare $U^*_\infty$ at $\theta = 0$ and $\theta = 180$ in figure 1), which produces azimuthally varying lift and drag forces. The net result is a larger upstream flow deceleration on the advancing side, i.e. the rotor appears less porous, leading to a larger wake deficit than on the retreating side. Generally, more porous bluff bodies produce weaker vortex shedding (Castro 1971; Steiros & Hultmark 2018; Steiros et al. 2020). However, the region of high wake deficit does not remain in the same position as the wake progresses downstream, but rather translates away from the centreline towards the advancing side (in the $+\hat {y}$ direction here). The second cause of wake asymmetry is a net cross-stream $-\hat {y}$ force on the blades (figure 6(d), dashed line) resulting in a net $+\hat {y}$ acceleration of the flow (figure 6e); here this cross-stream force is measured, but it may also be estimated as in Ayati et al. (2019). This net forcing is explained as follows: power measurements on a single-bladed cross-flow turbine (Strom et al. 2017) indicate that the majority, if not all of the power, is produced on the upstream side of the rotor, with peak power production centred at approximately $\theta = 90^\circ$ (rather than exclusively on the advancing side, as suggested by Hohman et al. 2018). This is illustrated by the location of the largest red tangential arrow in figures 6(a) and 12(a). As a consequence of this application of force, the fluid must experience a force in the opposite direction, specifically, in the $+\hat {y}$ direction. This effect is analogous to the angular velocity induced during axial-flow turbine operation, where the induced flow is opposite the direction of turbine rotation (Burton et al. 2001). In the case of the cross-flow turbine, the flow velocity induced in the cross-stream direction is advected downstream through the rotor to the wake, as depicted in figure 6(b). Strong evidence of this cross-stream velocity is seen in figure 5(c), though the action of blade tip vortices could also induce flow in this direction (Battisti et al. 2011).

Figure 6. (a) Measured streamwise, cross-stream and resulting tangential force vectors on a single-bladed turbine. Measurement methods and a demonstration of the validity of using single-bladed turbine measurements as a proxy for the force on one blade of a two-bladed turbine are given in Strom et al. (2017). (b) Average streamwise force on the blade as a function of cross-stream blade position ($0^\circ \leq \theta \leq 180^\circ$). Forcing on the advancing and retreating sides is compared, showing a larger streamwise forcing on the advancing side. (c) A cartoon of the effect on the wake. Flow is decelerated more heavily on the advancing side due to larger streamwise forcing, resulting in a larger wake deficit. (d) Cross-stream force on the blade as a function of azimuthal angle (solid) and average value (dashed). The average force is downward towards the retreating side. (e) An illustration of the resulting flow acceleration and effect on the wake: convection towards the advancing side, resulting in wake skew.

Returning to the properties of the mean wake, it is curious to note significant vertical (axial) velocities present in figure 5(d). Because we are sampling on the mid-plane and the turbine rotor is symmetric about this plane, one would expect the wake to reflect this vertical symmetry, resulting in no out-of-plane velocities. However, Peng et al. (2016) and Rolin & Porté-Agel (2015) both observe similar asymmetries. Interactions with the free-surface or the flume floor boundary layer could be mechanisms responsible for this phenomenon, although the former is unlikely as mid-plane vertical flows have been observed in wind-tunnel measurements. The stability of coherent wake structures may play some role in this asymmetry, as described later.

4. Results and analysis: periodic structures

Flows with natural or forced periodicity, such as the wake of a cross-flow turbine, contain turbulent fluctuations that are semi-regular in space or time, and thus differ from the stochastic fluctuations that occur further down the turbulent cascade. It is then useful to analyse flows with periodic, organized content in terms of the triple decomposition of Hussain & Reynolds (1970)

(4.1)\begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = \bar{\boldsymbol{u}}(\boldsymbol{x}) + \tilde{\boldsymbol{u}}(\boldsymbol{x},\phi(t)) + \boldsymbol{u}'(\boldsymbol{x},t), \end{equation}

where the total flow, $\boldsymbol {u}$, is the superposition of a time-averaged flow, $\bar {\boldsymbol {u}}$, the periodic flow, $\tilde {\boldsymbol {u}}$, parameterized by the phase $\phi (t)$ and incoherent fluctuations, $\boldsymbol {u}'$. While $\bar {\boldsymbol {u}}$ is calculated through simple time averaging, there are multiple approaches for separating the periodic and turbulent components. For flows where the forcing mechanism or flow periodicity are measured simultaneously with the velocity field, $\tilde {\boldsymbol {u}}$ can be calculated by phase averaging (i.e. by computing the ensemble mean of measurements occurring at the same phase of the forcing oscillator). This was the approach taken by the pioneers of the triple decomposition (Hussain & Reynolds 1970), and is commonly employed, for example in the wake of an axial-flow wind turbine by Eriksen & Krogstad (2017). Drawbacks to this approach include the necessity of measuring the phase of the forcing oscillator, as well as the introduction of statistical uncertainty in the case that measurements are not locked to the forcing oscillator (Cantwell & Coles 1983). This is the case for our measurements because PIV data collection was free running and not locked to the turbine blade position; however, PIV triggering was time synchronized with measurements of turbine performance and blade position. A potential solution for the uncertainty in free-running measurements is to use a weighted average based on the phase offset for a given measurement from the phase in question.

Alternatively, Fourier averaging, where $\tilde {\boldsymbol {u}}$ is estimated using a truncated Fourier series (Sonnenberger, Graichen & Erk 2000), eliminates the need to measure the forcing signal simultaneously with flow measurements and removes error associated with phase uncertainty. However, since the base forcing frequency must be known or assumed, periodic flow structures due to phenomena other than the primary forcing mechanism may not be included in $\tilde {\boldsymbol {u}}$.

The desire to automatically extract and rank the importance of spatially coherent and temporally periodic flow phenomena at multiple scales, without a priori knowledge of the frequencies of interest, has inspired a number of methods. Instead of an oscillatory component composed of a single base frequency, these methods yield a triple decomposition of the form

(4.2)\begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = \bar{\boldsymbol{u}}(\boldsymbol{x}) + \sum_{n=1}^R \tilde{\boldsymbol{u}}_n(\boldsymbol{x},\phi_n(t)) + \boldsymbol{u}'(\boldsymbol{x},t), \end{equation}

where $R$ is the number of oscillatory modes used in the reconstruction and $n$ is the mode number. Because the frequency, amplitude and phase of oscillations are determined directly from velocity time series, this potentially reduces errors inherent to conditional averaging of free-running data and discrete Fourier transform (DFT) methods.

The DMD (Rowley et al. 2009; Schmid 2010; Tu et al. 2014; Kutz et al. 2016) provides a scalable and data-driven approach to extract the periodic component $up_n$ in (4.2). DMD is a combination of proper orthogonal decomposition (POD) (Berkooz, Holmes & Lumley 1993; Holmes et al. 2012; Taira et al. 2017, 2020) in space and the Fourier transform in time. Although POD is widely used for spatial mode extraction (Oberleithner et al. 2011; Edgington-Mitchell et al. 2014), including in the wake of an axial-flow turbine (Lignarolo et al. 2015; Premaratne, Wei & Hu 2016), the modes are known to mix frequency content. This is illustrated by the fact that the snapshot POD modes (Sirovich 1987; Aubry et al. 1988; Brunton & Kutz 2019) do not depend on the order of the flow data time series. Many instances of the failure of POD to extract dynamically important modes for multi-scale systems have been documented (Sayadi et al. 2012; Baj, Bruce & Buxton 2015; Taira et al. 2017). In contrast, DMD modes are a linear combination of POD modes, specifically designed to be coherent in space and have distinct oscillation frequencies, as well as growth or decay rates. The recursive DMD algorithm of Noack et al. (2016) combines favourable aspects of both approaches, namely low residual prediction error, pure frequency content, orthonormality of modes and the interpretation of mode amplitudes as energy content, making it a valuable technique for reduced-order modelling.

4.1. Approaches for triple decomposition

Here, we describe several triple decomposition methods. In the next section, we will compare the efficacy of these methods for detecting oscillatory structures and their performance based on error and energy capture.

4.1.1. Blade position conditional averages

We compute three conditional averages based on the turbine blade position at the time of PIV image capture. In these methods, PIV data are binned based on blade position. Subsequently, the median, mean or weighted mean flow field velocities are calculated.

For a single point in space, all $n$ measurements are collected for which the blade position, $\theta$, satisfies

(4.3)\begin{equation} | \theta - \theta_i | < \Delta \theta, \end{equation}

where $\Delta \theta$ is the half-bin width and $\theta _i$ denotes the $i$th bin centre. The bin mean is

(4.4)\begin{equation} \tilde{\boldsymbol{u}}(\boldsymbol{x},\theta_i) = \frac{1}{n} \sum_{j = 1}^n (\boldsymbol{u}(\boldsymbol{x},\theta_j)-\bar{\boldsymbol{u}}(\boldsymbol{x})), \end{equation}

which makes use the property that the mean of the stochastic component is zero. In an effort to reduce the sensitivity of this method to potential measurement outliers, the bin median is similarly computed. The bin-mean method also introduces error through gradients in the flow field over the range of blade positions over the bin width. One solution is to shrink the bin width, but the number of measurements $n$ occurring in the bin vary inversely with $\Delta \theta$. If $n$ is too small, the stochastic fluctuations have a non-zero mean. To reduce bin-width error while maintaining higher statistical certainty, a weighted average is computed, where the weight of each measurement varies inversely with its distance from the bin centre

(4.5)\begin{equation} \tilde{\boldsymbol{u}}(\boldsymbol{x},\theta_i) = \frac{\displaystyle \sum_{j = 1}^n (\boldsymbol{u}(\boldsymbol{x},\theta_j)-\bar{\boldsymbol{u}}(\boldsymbol{x}))|\theta_i+\Delta \theta - \theta_j|}{\displaystyle\sum_{j = 1}^n |\theta_i+\Delta \theta - \theta_j|}. \end{equation}

In these methods, each half-rotation of the rotor is assumed to be one period of flow oscillation due to the symmetry of the two-bladed rotor. Reconstruction error, when compared with the full flow field, was minimized with $\Delta \theta = 3^\circ$, or 30 bins per half-revolution, resulting in, on average, $n = 67$ flow snapshots per bin.

As none of these approaches is entirely satisfactory in handling the error between the actual blade position and the position of the bin centre introduced by free-running data acquisition, this motivates an exploration of alternative methods.

4.1.2. Fourier series reconstruction

A Fourier-series-based reconstruction of $M$ harmonics is given by

(4.6)\begin{equation} \tilde{\boldsymbol{u}}(\boldsymbol{x},\theta(t)) = \sum_{m = 1}^M A_m(\boldsymbol{x}) \sin[m \omega_b t + \phi_m(\boldsymbol{x})], \end{equation}

wherein a series of sinusoidal functions are fit to the data. The resulting function is used to reconstruct the data at the bin-centre blade positions, eliminating the inter-bin blade position error of the previous method.

In the case of this flow, selection of the base oscillation frequency, $\omega _b$, is simple given that the blade passing frequency is the primary driver of flow oscillations. The flow field is computed, similar to the conditional-average methods, by reconstruction at times $t_i$ that correspond to bin centres $\theta _i$. In practice, the coefficients and phase fields $A_m(\boldsymbol {x})$ and $\phi _m(\boldsymbol {x})$ are computed via the windowed fast Fourier transform (FFT), with the time series padded appropriately to ensure the FFT output includes all $r$ frequencies exactly equal to $m \omega _b$, removing potential frequency interpolation error. This is referred to in the following sections as the DFT method. We evaluate two versions of this method. First, for the ‘DFTc’ method, the base oscillation frequency (the blade passing) is calculated from the average location of the largest peak of the flow data spectra. Second, in the ‘DFTm’ method, the base frequency is calculated from the encoder data collected during turbine operation.

4.1.3. Multi-modal decomposition via optimized DMD

In many cases, the base oscillation frequency may be unknown, or the flow may exhibit features that oscillate at unrelated frequencies. In the case of a cross-flow turbine wake, the blade passing frequency may not be the only mechanism determining the time scale of periodic fluctuations. This motivates a generalization of the triple decomposition to (4.2) as introduced by Baj et al. (2015). Here, $\tilde {\boldsymbol {u}}$ is split into fluctuating components whose frequencies are not necessarily related, allowing this data-driven triple decomposition method to be used as an exploratory/diagnostic tool. In this work, we used DMD to identify the fluctuating components. A related method that could be used is spectral POD, an implementation of the original POD of Lumley (1967), with the ‘spectral POD’ terminology introduced by Picard & Delville (2000). A detailed account of the relationship between spectral POD and DMD is given by Towne, Schmidt & Colonius (2018).

The DMD was introduced by Schmid (2010) in the fluids community to identify spatiotemporal coherent structures from time-series data. In its simplest form, the DMD algorithm extracts the dominant eigenvalues and eigenvectors of the best-fit linear operator that approximately advances the measured state forward in time. The DMD algorithm starts with two snapshot matrices constructed of spatial and temporal flow components

(4.7a,b) \begin{align} \boldsymbol{X} = \begin{bmatrix} u(\boldsymbol{x}_1,t_1) & u(\boldsymbol{x}_1,t_2) & & u(\boldsymbol{x}_{1},t_{m-1}) \\ \vdots & \vdots & & \vdots \\ u(\boldsymbol{x}_n,t_1) & u(\boldsymbol{x}_n,t_2) & & u(\boldsymbol{x}_n,t_{m-1}) \\ v(\boldsymbol{x}_1,t_1) & v(\boldsymbol{x}_1,t_2) & & v(\boldsymbol{x}_{1},t_{m-1}) \\ \vdots & \vdots & {{\cdot}\mkern1mu{\cdot}\mkern1mu{\cdot}} & \vdots \\ v(\boldsymbol{x}_n,t_1) & v(\boldsymbol{x}_n,t_2) & & v(\boldsymbol{x}_n,t_{m-1}) \\ w(\boldsymbol{x}_1,t_1) & w(\boldsymbol{x}_1,t_2) & & w(\boldsymbol{x}_{1},t_{m-1}) \\ \vdots & \vdots & & \vdots \\ w(\boldsymbol{x}_n,t_1) & w(\boldsymbol{x}_n,t_2) & & w(\boldsymbol{x}_n,t_{m-1}) \\ \end{bmatrix},\quad \boldsymbol{X}' = \begin{bmatrix} u(\boldsymbol{x}_1,t_2) & u(\boldsymbol{x}_1,t_3) & & u(\boldsymbol{x}_{1},t_m) \\ \vdots & \vdots & & \vdots \\ u(\boldsymbol{x}_n,t_2) & u(\boldsymbol{x}_n,t_3) & & u(\boldsymbol{x}_n,t_{m}) \\ v(\boldsymbol{x}_1,t_2) & v(\boldsymbol{x}_1,t_3) & & v(\boldsymbol{x}_{1},t_m) \\ \vdots & \vdots & {{\cdot}\mkern1mu{\cdot}\mkern1mu{\cdot}} & \vdots \\ v(\boldsymbol{x}_n,t_2) & v(\boldsymbol{x}_n,t_3) & & v(\boldsymbol{x}_n,t_{m}) \\ w(\boldsymbol{x}_1,t_2) & w(\boldsymbol{x}_1,t_3) & & w(\boldsymbol{x}_{1},t_m) \\ \vdots & \vdots & & \vdots \\ w(\boldsymbol{x}_n,t_2) & w(\boldsymbol{x}_n,t_3) & & w(\boldsymbol{x}_n,t_{m}) \\ \end{bmatrix}. \end{align}

The best-fit linear operator that maps $\boldsymbol {X}$ into $\boldsymbol {X}'$ is given by $\boldsymbol {A}$, satisfying the approximate relationship

(4.8)\begin{equation} \boldsymbol{X}'\approx \boldsymbol{A}\boldsymbol{X}. \end{equation}

In practice, this matrix $\boldsymbol {A}$ may be approximated using the pseudo-inverse of $\boldsymbol {X}$, which is computed by taking the singular value decomposition $\boldsymbol {X}=\boldsymbol {U}\boldsymbol {\varSigma }\boldsymbol {V}^T$ and inverting each of the matrices $\boldsymbol {U}$, $\boldsymbol {\varSigma }$, and $\boldsymbol {V}^T$

(4.9)\begin{equation} \boldsymbol{A} = \boldsymbol{X}'\boldsymbol{V}\boldsymbol{\varSigma}^{{-}1}\boldsymbol{U}^T. \end{equation}

The matrix $\boldsymbol {\varSigma }$ is diagonal, and both $\boldsymbol {U}$ and $\boldsymbol {V}$ are unitary, so their transposes are their inverses. However, if the state $\boldsymbol {X}$ is a large discretized fluid velocity or vorticity field, the matrix $\boldsymbol {A}$ may be intractably large to represent, let alone to analyse. Instead, we compute the projection of $\boldsymbol {A}$ onto the leading POD modes, given by the first $r$ columns of $\boldsymbol {U}$, denoted by $\boldsymbol {U}_r$

(4.10)\begin{equation} \tilde{\boldsymbol{A}} = \boldsymbol{U}_r^T\boldsymbol{A}\boldsymbol{U}_r = \boldsymbol{U}_r^T\boldsymbol{X}'\boldsymbol{V}_r\boldsymbol{\varSigma}_r^{{-}1}. \end{equation}

The matrices $\boldsymbol {A}$ and $\tilde {\boldsymbol {A}}$ share the same eigenvalues, so it is possible to compute the spectrum of $\boldsymbol {A}$ by computing the eigendecomposition of $\tilde {\boldsymbol {A}}$

(4.11)\begin{equation} \tilde{\boldsymbol{A}}\boldsymbol{W} = \boldsymbol{W}\boldsymbol{\varLambda}, \end{equation}

where $\boldsymbol {W}$ contain the eigenvectors of $\tilde {\boldsymbol {A}}$ and $\boldsymbol {\varLambda }$ contains the eigenvalues. Finally, it is possible to compute the high-dimensional eigenvectors $\boldsymbol {\varPhi }$ of the matrix $\boldsymbol {A}$ (e.g. the DMD modes), from the low-dimensional eigenvectors $\boldsymbol {W}$ using the exact DMD algorithm of Tu et al. (2014)

(4.12)\begin{equation} \boldsymbol{\varPhi} = \boldsymbol{X}'\boldsymbol{V}_r\boldsymbol{\varSigma}_r^{{-}1}\boldsymbol{W}. \end{equation}

DMD has recently been connected to spectral POD (Towne et al. 2018), used to analyse a cross-flow turbine wake by Araya et al. (2017), and the resolvent operator (Sharma, Mezić & McKeon 2016). Another view is that DMD is an approximation of the Koopman operator, which is an infinite-dimensional linear operator that steps a system forward in time by operating on an infinite-dimensional Hilbert space of all scalar-valued functions of system measurements (Rowley et al. 2009; Mezić 2013; Kutz et al. 2016).

It is well known that the original DMD algorithm of Schmid (2010) is sensitive to noise (Bagheri 2014) and there are several recent approaches to de-bias the algorithm for noisy data (Dawson et al. 2016; Hemati et al. 2017; Askham & Kutz 2018). The optimized DMD (optDMD) algorithm of Askham & Kutz (2018) considers the evolution of all of the snapshots at once, instead of through a single iteration through the map $\boldsymbol {A}$, and provides an efficient way of solving a nonlinear least-squares regression problem using variable projection. This has the added benefit of allowing for an optimal DMD fit from data that are unevenly spaced in time. This is, in general, a non-convex procedure, although there are efficient algorithms to compute this optimization, and the results indicate considerable noise robustness over standard algorithms.

The optimized DMD method also provides a mechanism for constraining the eigenvalues of the returned modes, for example to keep them on the unit circle. This allows for solving of periodic-only optDMD modes, and can be used to restrict the oscillation frequencies. The data taken in these experiments consist of overlapping fields of view taken at separate times. When optDMD is performed on the entire dataset, the resulting modal oscillations are out of phase. The field-of-view overlap regions are used to correct the phase misalignment, resulting in full-field DMD modes. This method is likely useful for modal analysis in any experiment utilizing multiple overlapping measurement areas. Details of this method can be found in Nair et al. (2020).

The ability to restrict DMD eigenvalues to lie on the unit circle is critical for application to the triple decomposition, where such behaviour is inherent in the definition of the oscillatory term. The standard DMD algorithm could be used to determine oscillatory flow components by either selecting modes with imaginary-only eigenvalues, or by manually zeroing the real part of the eigenvalues. However, the original mode shapes returned by exact DMD are no longer guaranteed to best represent the data given the now altered eigenvalues. Optimal DMD circumvents this issue by iteratively optimizing the mode shapes given constraints on the eigenvalues.

4.2. Decomposition method comparison

For each of the algorithms above, the periodic component is extracted, reconstructed for the full length of the original dataset, and then added back to the mean flow ($\bar {\boldsymbol {u}} + \tilde {\boldsymbol {u}}$). The reconstruction is compared with the original flow in two ways. First, the average $L_2$ error between the reconstruction and original data is computed. A smaller $L_2$ error indicates that more of the oscillatory mode information is being captured. However, as the original data contain stochastic fluctuations, this $L_2$ error is never identically zero; note that the stochastic fluctuations are assumed to have no periodicity, and thus are contained only in the third component of the triple decomposition. Second, the total sum of the flow kinetic energy over space and time is computed for the original and reconstructed flow. The ratio of these energies is indicated by the vertical axis of figure 7(a), while the error is on the horizontal axis. The DMD-based method results in a reconstructed flow field with more energy explained and a lower error. Somewhat surprisingly, the DFTm method results in an order of magnitude higher error than DFTc as small errors in the encoder-measured frequency vs the true frequency results in large oscillation phase errors during reconstruction. This illustrates the importance of knowing or calculating the base frequency of interest exactly when using a DFT-based method. Of the averaging methods, a bin median, which is resilient to outliers, explains more of the flow kinetic energy than a bin mean or weighted bin mean.

Figure 7. (a) Kinetic energy content of the mean plus the reconstructed periodic flow normalized by the kinetic energy content of the full flow measurements vs the $L2$ error of the reconstruction vs the original flow. We expect the most effective triple decomposition method to minimize the error while maximizing the amount of energy capture (as indicated by the arrow). (b) Power spectra of the modes of the DFT and DMD methods. DMD indicates importance of low-frequency modes that may not be discovered by other methods. The first four DMD modes are labelled $M 1\to 4$.

The DMD method does not require a priori knowledge of the base frequency, and can be used to uncover flow phenomena oscillating at related or unrelated frequencies. As illustrated by figure 7(b), the first seven modes extracted by the optDMD algorithm contain the blade passing frequency and its first harmonic, followed by five lower-frequency modes. To view these modes in terms of the full wake, not just the separate fields of view collected, the phase of oscillation of each mode in each individual field of view was adjusted via numerical optimization to match the oscillation of neighbouring fields of view. This process is illustrated in figure 8 and details are given in Nair et al. (2020).

Figure 8. DMD mode phase correction process shown on the first DMD mode of the turbine wake data. Because data were collected at differing times, the phase of oscillation of the same mode in differing fields of view are not aligned. A numerical minimization of the error in field-of-view overlap regions is used to correct the phase. For this example, this is a five variable optimization problem (one field of view is the reference).

4.3. Wake DMD modes

The optDMD triple decomposition extracts and ranks oscillatory modes in terms of energy content. The first four modes are shown in figure 9 and their corresponding frequencies are identified in figure 7(b). Modes are computed using all three velocity components as measurement inputs. For compactness, the horizontal velocity magnitude

(4.13)\begin{equation} U_H = \sqrt{u^2 + v^2}, \end{equation}

is plotted in figure 9. As mentioned previously, modes one and two correspond to the blade pass frequency and its first harmonic. The first harmonic energy content decreases faster in the downstream direction than the fundamental frequency, perhaps due to the more rapid dissipation of the smaller-scale features. The structures responsible for the energetic oscillations at the blade passing frequency will be discussed in the following section. The frequency of the third mode is half the blade passing frequency. Close to the turbine, energy in this mode illustrates changes in flow due to small geometric differences in the rotor blades or their mounting angle. Together, modes three and four illustrate a phenomenon on the advancing side of the wake. Structures that initially occur at the blade passing frequency seem to be breaking down or combining into lower-frequency structures in a repeatable manner. This could be evidence of a transition toward the bluff-body far wake oscillation documented by Araya et al. (2017). However, the frequency of mode four is 1.18 Hz, while the predicted bluff-body frequency for a cylinder of the same diameter as the rotor is 0.8 Hz. It is possible that full transition to the bluff-body frequency has not yet occurred, and that measurements made further downstream would show lower dominant frequencies.

Figure 9. Modes extracted using the optDMD algorithm, with the phase of oscillations corrected. Modes are ranked by energy content and are identified by the M labels in figure 7(b).

5. Results and analysis: coherent structures

To provide a detailed description of the dynamics of coherent structures in the non-stochastic component of the wake ($\bar {\boldsymbol {u}} + \tilde {\boldsymbol {u}}$), we present the finite-time Lyapunov exponent (FTLE) fields (Haller 2002; Shadden, Lekien & Marsden 2005; Green, Rowley & Haller 2007; Farazmand & Haller 2012), in particular the Lagrangian coherent structure (LCS) ridges of the FTLE field. FTLE detects coherent structures via their boundaries by computing the maximum strain rate of a Lagrangian packet of fluid over a finite time period. Integration backward and forward in time yields attracting and repelling FTLE ridges, respectively, which enclose coherent structures. Unlike Eulerian methods, such as the Q-criterion (Hunt, Wray & Moin 1988), the FTLE field does not require a user-defined, subjective threshold to identify coherent structures. Additionally, as an integration-based method, FTLE analysis is more robust to noise than derivative-based methods (Green et al. 2007). Although other methods have been shown to identify LCS more rigorously (Haller & Beron-Vera 2013), the ridges of the FTLE field remain practically useful due to their ease of computation, mathematical simplicity and conceptual accessibility. Here, structures associated with shear-layer roll-up, blade-vortex shedding and bluff-body wake oscillation are identified, and the role of these structures in wake mixing and recovery is considered. In addition, vortex-core tracking is performed on the full time-resolved wake to determine vortex longevity and trajectory repeatability, which has implications for interactions with downstream turbines in an array.

5.1. Computing the FTLE

Lagrangian coherent structures are useful for identifying coherent regions of unsteady fluid flows that are segmented by time-varying separatrices, which are the unsteady analogues of stable and unstable invariant manifolds in dynamical systems (Haller 2002; Shadden et al. 2005). LCS are often computed as the second derivative ridges of the FTLE field (Shadden et al. 2005), which describes the maximum rate of stretching of a Lagrangian packet of fluid over a finite time period. The more recent method of Farazmand & Haller (2012) uses variational theory to compute LCS from FTLE fields.

The FTLE field is generally computed by integrating passive tracer particles along the flow of the velocity field $\boldsymbol {u}(\boldsymbol {x},t)$ for a time span $T$ as

(5.1)\begin{equation} \varPhi_{t_0}^T(\boldsymbol{x}(t_0)) = \boldsymbol{x}(t_0+T) = \boldsymbol{x}(t_0) + \int_{t_0}^{t_0 + T} \boldsymbol{u}(\boldsymbol{x}(\tau),\tau)\,{\rm d} \tau, \end{equation}

where $\varPhi _{t_0}^T$ is the flow map. Next, the flow map Jacobian, $D\varPhi _{t_0}^T$ is approximated via finite-difference derivatives with neighbouring points in the flow. In two dimensions, the flow map Jacobian at a point $\boldsymbol {x}_{i,j}$ is

(5.2) \begin{align} (D\varPhi_{t_0}^T)_{i,j} \approx \begin{bmatrix} \dfrac{x(t_0 + T)_{i+1,j}-x(t_0 + T)_{i-1,j}}{x(t_0)_{i+1,j}-x(t_0)_{i-1,j}} & \dfrac{x(t_0 + T)_{i,j+1}-x(t_0 + T)_{i,j-1}}{y(t_0)_{i,j+1}-y(t_0)_{i,j-1}}\\ \dfrac{y(t_0 + T)_{i+1,j}-y(t_0 + T)_{i-1,j}}{x(t_0)_{i+1,j}-x(t_0)_{i-1,j}} & \dfrac{y(t_0 + T)_{i,j+1}-y(t_0 + T)_{i,j-1}}{y(t_0)_{i,j+1}-y(t_0)_{i,j-1}} \end{bmatrix}. \end{align}

From a continuum mechanics standpoint, this is a numerical computation of the deformation gradient. The FTLE $\sigma$ is computed from the largest eigenvalue $\lambda _{max}$ of the Cauchy–Green deformation tensor $\boldsymbol {\varDelta }=(D\varPhi _{t_0}^T)^\intercal D\varPhi _{t_0}^T$, which is the maximum singular value of the flow map Jacobian

(5.3)\begin{equation} \sigma(\boldsymbol{x}_0,t_0,T) = \frac{1}{T}\ln(\sqrt{\lambda_{max}[\boldsymbol{\varDelta} (\boldsymbol{x}_0,t_0,T)]}). \end{equation}

The variable $\sigma$ is a scalar field that is typically computed on a discrete grid of particles, and for unsteady flows this field is recomputed at every time step, introducing redundant computations that may be eliminated (Brunton & Rowley 2010; Luchtenburg, Brunton & Rowley 2014). When $\sigma$ is large, then neighbouring particles undergo considerable stretching along the flow, while particles with small $\sigma$ will tend to remain in coherent patches with their neighbours. Repelling or attracting FTLE field structures may be computed by integrating particles either forward or backward time, respectively.

As the integration period, $T$, is increased, the ridges of the FTLE field become more defined, although their locations remain constant (Green et al. 2007). For studies with limited interrogation windows, such as this one, this can present a problem. Increasing integration time results in more finely resolved FTLE ridges, but increases the chance that a passive tracer will exit the flow field before the end of the integration. A partial solution to this issue is to subtract the global (time and space) mean velocity vector from the flow field, as this will not affect the mechanics of a Lagrangian fluid packet. However, with flows that exhibit local velocities differing significantly from the global mean, this strategy is only partially effective since the local velocities can still eject the passive tracers. FTLE computations presented here assign passive tracers that leave the domain their FTLE value at the time of exit, resulting in a decrease in FTLE ridge definition near some domain edges.

5.2. Wake coherent structures

In addition to analysing the frequencies of different components of the oscillating wake, it is illuminating to examine the time evolution of periodic coherent structures. In figure 10 we show the forward- and backward-time FTLE fields computed from the mean and the reconstructed periodic component (optDMD method). We observe the structures driving the flow oscillation that were evident in the first two DMD modes. On the advancing side of the wake ($y/D > 0$), there is a vortex street with vorticity opposite the direction of turbine rotation (figure 11). Roll up of the strong shear layer in the flow, as shown in figures 5(a) and 5(b), likely energizes and sustains these vortices. However, their initial source may be the starting vortex shed by the blade as it passes through $\theta = 0^\circ$ and lift production for the upstream side of the stroke commences (see upper half of figure 12a). Alternatively, the disturbance caused by the blade passage could be enough for the strong shear layer to roll up at a regular frequency, regardless of lift production.

Figure 10. Forward and backward FTLE fields computed on $\bar {\boldsymbol {u}} + \tilde {\boldsymbol {u}}$ (optDMD method). These fields represent areas of maximum stretch and convergence, respectively, and together outline the boundaries of coherent structures.

Figure 11. Second derivative ridges of the FTLE field superimposed on the out-of-plane flow vorticity (the curl of the horizontal velocity components).

Figure 12. (a) Measured power coefficient as a function of azimuthal blade position for a single blade in a cross-flow turbine. (b) Vorticity shed as a result of lift generated during the primary power production region of the blade stroke (red region in (a)).

On the retreating side of the wake ($y/D < 0$), a larger vortex structure is apparent. During the power-producing portion of the blade stroke (centred around $\theta = 90^\circ$ when the blade is farthest upstream) lift production requires the generation of counter-clockwise circulation around the blade. Because this circulation is not permanent, but is created every rotation, packets of opposite circulation must be shed into the wake (Battisti et al. 2011). A dynamic stall event could force a sudden shedding of this circulation, resulting in a much stronger structure than on the advancing side, where lift increases comparatively gradually, possibly dispersing the starting vortex (see figure 12). Detailed analysis of the actual formation of this vortex structure on the blade would require further measurements upstream. This structure dissipates more quickly than the vortex street on the advancing side, perhaps due to the lack of a strong energizing shear layer on this side of the wake. Alternatively, the vortex street may be the source of the weaker shear layer due to cross-stream mixing caused by this structure. Additionally, cross-stream mixing may contribute to the faster wake recovery on the retreating side, and thus to the overall asymmetry of the mean wake (see figure 5a).

The vortex structure shed on the retreating side is more complex than a single rotating packet of fluid. It appears that an intense core of vorticity is surrounded by a ring of vorticity of the same rotation direction. This unstable arrangement may be partially responsible for the rapid breakdown of this structure, although to see why it is able to survive at all we must examine the out-of-plane ($w$) velocity, shown in figure 13. The inner vortex core has intense axial (vertical) velocity in the negative direction, while the outer vorticity ring has vertical velocity towards the flume free surface. We see now that the out-of-plane velocity observed in the mean wake (figure 5d) is entirely due to this structure. It appears that this structure already contains significant axial flow when it is shed into the wake, so it seems unlikely that it is due purely to asymmetries in the free-stream velocity and their influence on shed tip vortices. Axial flow in dynamic-stall-related structures has been reported in insect flight (Birch & Dickinson 2001) and delta wing aircraft (Wu, Vakili & Wu 1991). However, the source of the pressure gradient that can drive vertical flow within a rotor that is symmetric about the mid-plane remains unknown.

Figure 13. Out-of-plane (vertical) velocity, mean and periodic component ($\bar {w} + \tilde {w}$).

The predictability of the trajectories of wake coherent structures is of interest for optimizing the performance of closely spaced downstream turbines in an array. Consistent trajectories may make it easier for a downstream turbine to harness or avoid coherent structures (i.e. operation can be coordinated based on knowledge of upstream turbine phase). In figures 14(a) and 14(b), the core of the retreating side vortex is tracked. All tracks are shown in figures 14(c) and in 14(d), the probability density of the cross-stream location of the vortex core is plotted as a function of streamwise vortex-core position. The position of the vortex core increases in uncertainty rapidly with $x/D$. This indicates that the longevity of vortices may be greater than that predicted by just the periodic component of the flow, on which both the optDMD modes and, therefore, the subsequent FTLE are based. However, regardless of longevity, the rapid increase in the uncertainty of trajectories means that the structure could not be intercepted reliability by a downstream turbine blade for $x/D >1$.

Figure 14. Retreating side vortex-core tracking. (a,b) Example tracks and the corresponding vorticity fields. (c) All 50 tracks for the sampling period. (d) Probability distribution of track $y/D$ location.

6. Discussion and conclusions

This work has presented three detailed analyses of a cross-flow turbine wake based on planar PIV measurements of three-component velocities: (i) analysis of time average and relation to forces acting on the turbine; (ii) demonstration of optimized DMD to identify energetically important modes that cannot be discovered by other methods; and (iii) quantification and tracking the form and trajectory of coherent structures shed into the near wake to identify the region over which their propagation is deterministic, which is of relevance to array control. First, we show that the mean wake structure is similar to prior investigations, despite the relatively higher $c/R$ for this turbine. Using phase-resolved performance measurements, we hypothesize the observed wake skew is simply a consequence of momentum conservation for the torque-producing tangential force acting on the blade, clarifying the conflicting and inconsistent explanations in prior work. Second, we show that the optDMD algorithm can identify the periodic flow in a triple decomposition with lower error than conditional averaging or DFT-based methods, explains more of the energy present, and discovers oscillating structures at unrelated frequencies. As for Araya et al. (2017), we observe indications of a transition to bluff-body shedding in the far wake, particularly on the advancing side of the rotor. Third, building on the identified periodic structures, we present the first detailed description of Lagrangian coherent structures in a cross-flow turbine wake. Vortex streets on the advancing and retreating sides are observed and their formation mechanisms hypothesized. We observe remarkably high axial flow in the core of the vortices associated with lift production on the retreating side of the wake and note that this axial flow appears to originate within the confines of the rotor. The presence of a mean vertical flow through the mid-plane of a spanwise symmetric rotor motivates future investigation of its origins and connecting the near-wake structures observed here to the fluid–structure interactions at the blade.

As described in § 1, one of the factors driving modern interest in cross-flow turbine fluid dynamics is the potential for superior performance in arrays, particularly arrays of closely spaced turbines where interaction with mean and periodic coherent structures in one turbine's wake can be exploited by another. As shown here, the mean and periodic wake structure is temporally and spatially rich, with varying implications for turbine arrays. For example, coherent structures shed on the retreating side of the wake likely dissipate too quickly to interact with a turbine blade more than one diameter downstream. The coherent structures on the advancing side wake are more persistent and transition to lower-frequency structures. As these structures are similar to the rotor size and located in the slowest mean flow, they would likely be detrimental to downstream turbine performance at greater distances. We do note, however, while the coherent structures are visually striking (see figure 10), the energy contained in these structures is quite low (cf. magnitude of DMD mode 4 in figure 9). An interesting avenue of future work would be to explore other coherent structure extraction techniques, such as the acceleration feature points method (Kasten et al. 2016), which is well suited for vortex and saddle point extraction and the automated tracking of trajectories. Similarly, the FTLE results are based on DMD-aligned flow fields, which was necessary to analyse the full field of view across the six smaller PIV fields of view. However, it is an interesting avenue of future work to explore how the FTLE computation is affected by this DMD filtering; this would require a benchmark velocity field with a large field of view, where it is possible to compute FLTE before and after DMD filtering. Finally, we note that this dataset is at a Reynolds number that is lower than in many commercial applications such that an increase in Reynolds number may increase the dissipation rate of coherent structures, as well as decrease consistency in their trajectory and strength (Rocha, Bernardo & Franck 2018).