2.1. Wind tunnel and wing

The experiments were conducted in the large-scale closed-loop wind tunnel at the Department of Energy and Process Engineering of the Norwegian University of Science and Technology (NTNU). The wind tunnel has a test-section with dimensions of $1.8\,\mathrm{m}\,\times 2.71\,\mathrm{m}\,\times 11.1\,\mathrm{m}$ in height, width and length, respectively. It houses a six-component force balance that can be rotated to vary the wing model’s angle of attack and to obtain an accurate measure of the forces acting on it. Further details on the force balance measurements are given in § 2.3. The wind tunnel can also be fitted with an active turbulence grid $0.8\,\mathrm{m}$ from the start of the test section (after the tunnel’s convergent section), which enables the generation of highly turbulent flows. The particular active grid used in this facility is described in detail by Kildal et al. (Reference Kildal, Li, Hearst, Petersen and Øiseth2023). It spans the entirety of the test section’s cross-sectional area and consists of 90 shaft assemblies, with each assembly being individually operated by a dedicated integrated stepper motor. The mesh length of the grid is $M_g = 0.1$ m. The grid was operated with different ‘fully random’ (Hearst & Lavoie Reference Hearst and Lavoie2015) sequences whereby the grid is sent top-hat distributions of the rotational velocity, rotation period and acceleration for each shaft individually. In this study, only the rotational velocity was varied between cases as this parameter was shown to have the greatest impact on the produced turbulence levels (Hearst & Lavoie Reference Hearst and Lavoie2015).

The employed finite-wing model was a rectangular, square-tipped half-wing with a NACA 4412 profile. It has a chord length of $c = 0.25\,\mathrm{m}$ , a maximum thickness of $\tau = 0.03\,\mathrm{m}$ and spans $b=0.9\,\mathrm{m}$ from the floor to the centre of the test-section. This profile was chosen because it is a commonly used reference geometry for which comparative measurements exist in the literature (Kay et al. Reference Kay, Richards and Sharma2020; Cruz Marquez et al. Reference Cruz Marquez, Monnier, Tanguy, Couliou, Brion and Dupont2021b ). The airfoil was machine-cut from a single block of polyurethane (Ebaboard 0600) which is an isotropic material with well-known characteristics and properties (flexural modulus of $900 \pm 200$ MPa). The airfoil was deliberately made of an isotropic material so that there is no preference for a particular loading. The airfoil was mounted vertically in the wind tunnel test-section and centred on the force balance. The leading-edge was $3\,\mathrm{m}$ ( $= 30M_g = 12c$ ) downstream of the active grid as shown in the schematic in figure 1. The wing was rigidly mounted and clamped to the force balance at the root such that the rotation axis of the force balance passed through its midchord. This effectively made a pin joint, whereby the wing span and tip were free while the root was fixed. The mounting platform of the force balance was below the wind tunnel floor. Thus, only the NACA 4412 section was exposed to the flow and not the root or the clamping apparatus. Additionally, a $4\,\mathrm{mm}$ gap was left between the wind tunnel floor and the wing model so as to avoid contact, which might influence the force measurements.

The coordinate system used in this study, $ [ x,y,z ]$ , is as shown in the schematic in figure 1. The $x$ -axis is along the streamwise direction, the $y$ -axis along the transverse direction (width of the wind tunnel) and the $z$ -axis along the wing’s span (height of the wind tunnel). The corresponding velocity components are represented by $ [ u,v,w ]$ . The origin is at the wing’s tip at the trailing edge when its angle of attack is $\alpha = 5^{\circ }$ , given that most of the results presented and discussed herein correspond to this wing configuration. A second coordinate system based on the wingtip vortex’ centre will also be used herein, which is represented by $ [\hat {y},\hat {z} ]$ . In the rest of the document, $\overline {\boldsymbol{\cdot}s }$ represents time-averaged quantities, $\boldsymbol{\cdot}s '$ represents standard deviation, $\hat {\boldsymbol{\cdot}s }$ represents time-averaged quantities in the ( $\hat {y},\hat {z}$ ) coordinate system, which are sometimes referred to as conditionally averaged quantities, and $\langle \boldsymbol{\cdot}s \rangle$ represents azimuthal averaging.

Figure 1. Schematics of the experimental set-up (not to scale), showing the wind tunnel, active grid, flexible wing model with speckle pattern, turntable, force balance and the digital image correlation (DIC) and stereoscopic particle image velocimetry (SPIV) systems (including cameras, laser sheet and light-emitting-diode (LED) illumination). The coordinate system is also indicated. (a) Side view of the wind tunnel’s central plane (grey cameras were mounted on the tunnel’s side); (b) top view (grey SPIV camera located beneath the tunnel).

2.2. Test cases and incoming flow conditions

The experiments were conducted at a mean freestream velocity of $U_{\infty } \approx 7.8\,\mathrm{m\,s}^{-1}$ , corresponding to a chord-based Reynolds number ${\textit{Re}}_{c} = U_{\infty } c / \nu = 1.4 \times 10^{5}$ . This is within the range of Reynolds numbers investigated in previous studies (see table 1) and was selected because it was the minimum velocity whereby Reynolds number independence was observed.

Four turbulent flow conditions were generated and investigated for this study: REF, A, B and C. The first test condition, REF, was a low turbulence reference case without the active grid. For case A, the active grid was in place but kept static and fully open, producing a turbulent flow with low intensity. For the last two cases, the active grid was actuated in two modes to generate approximately homogeneous, isotropic turbulence. Case B involved spinning all the vertical axes of the active grid at a frequency of $7\pm 2$ Hz, while in the last case (C), both the horizontal and vertical axes were spun at a frequency of $5\pm 2$ Hz.

A preliminary study using HWA was carried out to characterise the incoming flow conditions corresponding to the four test cases and to ascertain flow homogeneity in the wind tunnel test section without the wing. A Dantec 55P11 single hot-wire probe, controlled by a Dantec StreamLine Pro constant temperature anemometer, was employed. The probe was fixed in close proximity to a Pitot-static tube during all the tests. This system was fixed to a four degree-of-freedom wind tunnel traverse in order to measure the incoming flow at different locations in the wind tunnel. A thermocouple was used to measure the flow temperature. The hot-wire signals were sampled at a nominal acquisition rate of $75\,\mathrm{kHz}$ , with an analogue low-pass cutoff filter set to $30\,\mathrm{kHz}$ . Before and after each test case, the hot-wire probe was calibrated against the Pitot-static tube at the wind tunnel’s centreline, $80M_g$ downstream of the active grid. The calibration was performed for 11 velocities (from $2.3$ to $15.5\,\mathrm{m\,s}^{-1}$ ). The temperature correction methodology of Hultmark & Smits (Reference Hultmark and Smits2010) was employed. During postprocessing, a seventh-order low-pass Butterworth filter was applied to the time series at $1.1 f_{\eta }$ , where $f_{\eta }$ is the Kolmogorov advection frequency.

To assess the flow characteristics in the wind tunnel for the four test cases, the base flows without the wing were evaluated along the centreline of the wind tunnel from $0.8c$ upstream of the leading edge of the wing to $2.2c$ downstream of the trailing edge in steps of $2M_g$ , as well as in a plane consisting of $3 \times 3$ points situated at the leading edge. For these tests, hot-wire samples for cases A, B and C were acquired for $120\,\mathrm{s}$ ; data were only acquired for $30\,\mathrm{s}$ for the REF case as this flow was nearly laminar. Specifically at $30M_g$ downstream of the active grid (the wing’s leading edge), data were acquired for $600\,\mathrm{s}$ for cases A, B and C, and $120\,\mathrm{s}$ for case REF, respectively, to converge the turbulent spectra. The mean velocity remained relatively constant for all test conditions, with standard deviations from the average value of the mean velocities along the centreline of 0.8 % for the worst case. The streamwise scan showed that the turbulence intensity ( $u'/\overline {u}$ ) remained around $0.2\,\%$ for REF, decayed from $3\,\%$ to $2.5$ % for A, from $8.5\,\%$ to $7.8\,\%$ for B, and from $13\,\%$ to $11.5\,\%$ for C. Such a result is expected as grid-generated turbulence decays in space in a wind tunnel (see, for example, Comte-Bellot & Corrsin Reference Comte-Bellot and Corrsin1966) and is reported here only for completeness.

A measure of the homogeneity of the flow along the transverse direction ( $y$ – $z$ plane) was obtained by conducting planar scans at $30M_g$ downstream of the active grid (the leading edge of the wing) were performed over a $4M_g \times 4M_g$ grid oriented perpendicular to the flow. The analysis array consisted of nine points, with the central point being the second point used in the centreline homogeneity analysis reported previously. The homogeneity in the mean and fluctuating velocity components of the flow over the $4M_g \times 4M_g$ array is shown in figure 2. The cases generally show homogeneity of the mean velocity within $\pm 2\,\%$ of the centreline velocity and reasonable uniformity of the turbulence intensity.

Figure 2. Homogeneity of the flow in a $4M_g \times 4M_g$ array $30M_g$ downstream of the active grid. The values for $\overline {u}$ (a) are normalised by the respective centreline velocity $\overline {u}_{(0,0)}$ . The values for $u'/\overline {u}$ (b) are the turbulence intensities as a percentage.

The turbulence properties of the incoming flow measured at a distance of $30M_g$ downstream of the inlet on the centreline for the four test cases are summarised in table 2. The table includes values for the turbulence intensity $u'/\overline {u}$ , the integral length scale $L_{uu}$ obtained from integrating the autocorrelation function of the velocity fluctuations to the first zero-crossing, the Taylor-microscale-based Reynolds number ${\textit{Re}}_{\lambda }{ = \lambda _T u'/\nu }$ and the global isotropy ratio. The square of the Taylor microscale is given by $\lambda _T^2 = u'^2/\langle (\partial u' / \partial x)^2\rangle$ and was calculated by estimating the gradients using Taylor’s hypothesis and a seventh-order centred-difference scheme as described in Hearst et al. (Reference Hearst, Buxton, Ganapathisubramani and Lavoie2012). Using the method described by Benedict & Gould (Reference Benedict and Gould1996), the maximum random error in the turbulence intensity measurement was determined to be approximately $1.5\,\%$ . The global isotropy ratio was measured by spatially averaging the isotropy ratio over the SPIV field-of-view (FOV) (described in § 2.4) acquired at $37.5M_g$ without the wing model present. Here, the transverse velocity components are averaged when put in the global isotropy ratio to create a single parameter. The isotropy ratio approaches unity as the turbulence level of the active grid is increased; note REF is not acquired with the active grid and thus does not follow the aforementioned trend and instead illustrates the background isotropy in the empty facility.

Table 2. Summary of the turbulence properties in the wind tunnel without the wing. Statistics are for $30M_g$ downstream of the active grid along the centreline ( $(y,z)=(0,0)$ in figure 2, wing’s leading edge), except for the isotropy ratios which were computed from PIV in the empty tunnel at the same location as the SPIV measurements for the wingtip vortex analysis, i.e. $37.5M_g$ from the grid or $2c$ downstream of the airfoil trailing edge.

The cases reported in table 2 illustrate that the FST intensity and ${\textit{Re}}_\lambda$ increase from case REF to case C; however, so does the integral length scale. As such, when the FST is described as ‘increasing,’ this should be understood to mean that the scale of the turbulence is increasing, i.e. it is both more energetic and contains energy at larger length scales.

The turbulent velocity spectra obtained from the hot-wire measurements for all cases are shown in figure 3. When viscous unit normalisation is used (figure 3 b), all turbulent spectra collapse at the small scales and differ only at the large scales, as one expects. A $\kappa ^{-5/3}$ line (black dashed) is added for reference. The absence of low-frequency peaks in the spectra indicates that the active grid forcing does not preferentially energise a particular frequency or set of frequencies, which is a desired outcome. As the turbulence intensity increases, the inertial range grows as expected from previous active grid measurements (Larssen & Devenport Reference Larssen and Devenport2011; Hearst & Lavoie Reference Hearst and Lavoie2015) and turbulence theory for increasing ${\textit{Re}}_\lambda$ .

Figure 3. Velocity spectra for all cases plotted in (a) dimensional frequency-space and (b) wavenumber-space normalised with viscous units. The REF case is only plotted in (a) because its turbulence level is too low to be meaningfully represented by the normalisation in (b).

2.3. Force measurements

The different forces acting on the wing model at various test conditions were measured using the six-component force balance that is housed in the wind tunnel. Before conducting force measurements, load cell calibration was performed using known standard reference weights. Offset measurements were taken before and after each set of measurements to monitor and account for any sensor drift. Linear interpolation was carried out to compensate for any drift. Force measurements were acquired for angles of attack in the range $-5^{\circ } \leqslant \alpha \leqslant +26^{\circ }$ . Data acquisition was performed at a sampling frequency of $2\,\mathrm{kHz}$ for each angle of attack. For cases REF and A, one $240$ s sample was acquired for each $\alpha$ . In order to obtain converged results for the two most turbulent cases, B and C, $540$ s of data were acquired at each $\alpha$ ; these data were acquired as three independent 180 s samples. By propagating the following sources of uncertainty – an assumed angle-of-attack uncertainty of $0.5^{\circ }$ and the residual square errors from the calibration – the maximum uncertainty in the measured lift force is estimated to be approximately $0.67\,\mathrm{N}$ , which corresponds to a relative uncertainty of approximately $5.7\,\%-9.5\,\%$ .

2.4. Stereoscopic particle image velocimetry measurements

Cross-stream SPIV was employed to measure the three components of the instantaneous velocity field in the wake of the wing. The FOV is in the $yz$ -plane. During these measurements, the wing was maintained at an angle of attack of $\alpha = 5^{\circ }$ . A preliminary smoke visualisation exercise was carried out at the four incoming turbulence conditions to estimate the tip vortex’ location, in order to place the SPIV system. Figures 4(a) and 4(b) show the evolution of the wingtip vortex downstream of the trailing edge for case C (the most turbulent case) at low velocity; low velocity was used to facilitate the qualitative visualisation. Based on this, the streamwise station $2c$ downstream of the wing’s trailing edge was chosen for the SPIV measurements as the vortex tended to remain within the plausible particle image velocimetry (PIV) measurement domain and there are reference measurements from previous experimental studies (Serrano-Aguilera et al. Reference Serrano-Aguilera, García-Ortiz, Gallardo-Claros, Parras and del Pino2016; Bölle et al. Reference Bölle, Brion, Couliou and Molton2023). The imaging set-up is shown in figure 1.

Figure 4. Smoke visualisation of the wingtip vortex (a) near the trailing edge, and (b) downstream up to approximately $30\,c$ . The flow case is C in these images and the velocity was reduced to decrease the dispersion of the smoke plume.

The air flow was seeded with di-ethyl-hexyl-sebacate particles of approximately $1\,\unicode{x03BC}\mathrm{m}$ diameter. The particles were illuminated with a Litron LDY303HE dual-pulse Nd:YLF laser emitting green light at a wavelength of $527\,\mathrm{nm}$ . The laser head was placed outside the wind tunnel. Using an appropriate set of laser optics, the beam was shaped into a ${\approx}1\,\mathrm{mm}$ -thick sheet, which entered the test section from the side through the available optical axis.

Imaging was carried out using two high-speed Phantom v2012 $1\,\mathrm{MP}$ cameras featuring a $1280 \times 800\,\mathrm{px}^2$ complementary metal-oxide-semiconductor (CMOS) sensor ( $28\,\unicode{x03BC} \mathrm{m}$ pixel pitch, $12\,\mathrm{bits}$ of digital resolution). These cameras were equipped with Nikon Nikkor $200$ mm macrolenses operated at $f_{\#}=5.6$ . Green filters were attached to the lenses to reduce the effect of other wavelengths (particularly from the LEDs used for DIC, see § 2.5) and improve the signal-to-noise ratio of the acquired SPIV images. The first camera was located on the side of the wind tunnel and directed towards the upstream side of the measured FOV, while the second camera was placed along the centreline ( $y = 0$ ) beneath the wind tunnel, downstream of the wing. In order to coincide with the measurement plane and to ensure that both cameras capture sharp and in-focus images of the particles in the same plane, Scheimpflug adapters were used to tilt the focal plane of each camera.

The laser, the two SPIV cameras and the two cameras used for DIC (described in § 2.5) were synchronised using a LaVision PTU X controller and the Davis 10.2 software suite. Images were captured in double-frame mode, with a pulse separation of $95\,\unicode{x03BC} \mathrm{s}$ between an image pair. For each test case, $10\,000$ image pairs were captured at a rate of $200\,\mathrm{Hz}$ , leading to a total sampling time of $50\,\mathrm{s}$ . Image preprocessing, calibration and vector computations were also carried out using Davis 10.2. Calibration was done using a LaVision Type 20 two-level calibration plate. Cross-correlation for vector field computation was carried out using multiple passes with a final interrogation window of $32\,\mathrm{px}\times 32\,\mathrm{px}$ and a relative overlap of $50\,\%$ , resulting in a vector spacing of approximately $3\,\mathrm{mm}$ in both the $y$ and $z$ directions. The final field was cropped to $54 \times 54$ vectors ( $160\,\mathrm{mm} \times 160\,\mathrm{mm}$ ) near the centre of the FOV where the vector quality was high.

Measurement uncertainties are quantified for the REF case. The estimated errors in the instantaneous in-plane velocity components ( $v,w$ ) are approximately $0.2\,\mathrm{m\, s}^{-1}$ , while the out-of-plane velocity component ( $u$ ) has an estimated error of $0.4\,\mathrm{m\, s}^{-1}$ . Given a freestream velocity of $7.8\,\mathrm{m\, s}^{-1}$ , these correspond to relative uncertainties of approximately $2.6\,\%$ and $5.1\,\%$ , respectively. Using the linear error propagation technique summarised by Sciacchitano & Wieneke (Reference Sciacchitano and Wieneke2016), the uncertainty of the mean velocity for the out-of-plane component (which has the highest uncertainty) was $\pm 0.4,\,0.9,\,0.9$ and $1.2\,\%$ for cases REF, A, B and C, respectively.

2.5. Digital image correlation measurements

Digital image correlation is an optical method that enables non-contact measurement of surface deformation and displacement. To this effect, a speckle pattern was applied on the wing, which can be seen in figure 4(a). The speckle pattern consisted of randomly arranged points, each approximately $0.5$ mm in diameter. Images of this pattern were acquired with no load to get the reference position. Then images were captured simultaneously with the PIV measurements to get a time series of the displacement and deformation. The experimental set-up used for the DIC is shown in figure 1. Imaging was carried out using two Photron SA 1.1 cameras featuring a $1024\,\mathrm{px} \times 1024\,\mathrm{px}$ monochrome CMOS sensor ( $20\,\unicode{x03BC}\mathrm{m}$ pixel pitch, $12\,\mathrm{bits}$ of digital resolution). They were positioned on the side of the wind tunnel and imaged the pressure side of the airfoil. The cameras were equipped with Sigma $50{-}100\,\mathrm{mm}$ lenses, which were adjusted to a numerical aperture of $f_{\#} = 1.8$ , and had a focal length of $100\,\mathrm{mm}$ . This resulted in a scaling of ${\approx}0.3\,\mathrm{mm\,px}^{-1}$ . The selected FOV was a rectangular area of $53.7\,\mathrm{mm}\times 218.4\,\mathrm{mm}$ in height and width, respectively, located at the tip of the wing. Note that the full sensor of the camera was not used. To achieve uniform illumination across the wing’s surface-of-interest, the set-up incorporated two GS Vitec LED light sources positioned between the two DIC cameras as shown in figure 1. Polarisers were used on both the LEDs and cameras to reduce glare. Finally, given the simultaneous operation of the SPIV and DIC systems, a $200\,\unicode{x03BC}\mathrm{s}$ delay was added to the DIC system to minimise the green laser’s interference with the DIC measurements.

Stereoscopic image calibration for DIC was done using a calibration cylinder with a reference pattern. A global DIC approach was employed to extract displacement information from the recorded images. A mesh with Q4 elements and an element size of 15 data points was used during processing. Bicubic spline interpolation and grey level normalisation were also employed. Images acquired during the REF case where the wing’s deformations were negligible (discussed in § 3.2) were used to quantify the noise level in the deformation measurements. The noise level was taken as the standard deviation of the deflections ( $y$ -direction deformations) over the entire FOV and averaged over 10 images. It was found to be $0.005\,\mathrm{mm}$ .

Figure 5. Impact of FST on the time-averaged (a) lift and (b) drag of the NACA 4412 flexible finite wing. The standard deviation of the (c) lift and (d) drag measurements that represent the fluctuations in the measured forces is also presented.

3.1. Forces on the wing

To understand how the finite-wing as a whole behaves in the produced turbulent flows, force balance measurements as introduced in § 2.3 were performed for $-5^{\circ } \leqslant \alpha \leqslant +26^{\circ }$ . As is typical of airfoil measurements, hysteresis was observed when comparing results for increasing and decreasing $\alpha$ over the same range. As such, all results presented herein are for increasing $\alpha$ . Figure 5(a) shows the variation of the time-averaged three-dimensional (3-D) lift coefficient with angle of attack for the four test cases. The 3-D lift coefficient is expressed as $C_{L} = L / ((1 /2) \rho U_{\infty }^{2} A )$ , where $L$ is the lift force measured by the force balance, $\rho$ is the air density and $A$ is the reference surface area of the wing. In the present results, slopes of the lift curves are approximately collapsed within the uncertainty of the measurement for the linear range. This suggests that FST has little influence on the mean aerodynamic behaviour of the finite wing in this region and is in line with the observations of Wang et al. (Reference Wang, Zhou, Mahbub Alam and Yang2014). A slight increase in the maximum lift coefficient is observed when increasing turbulence intensity. This is consistent with the works of Li & Hearst (Reference Li and Hearst2021) and Thompson et al. (Reference Thompson, Biler, Symon and Ganapathisubramani2023) but in opposition with the findings of Kay et al. (Reference Kay, Richards and Sharma2020), who observed a decrease in the maximum 2-D lift coefficient of a NACA 4412 airfoil when increasing the FST intensity from $1.3\,\%$ to $15\,\%$ .

From $\alpha = 15^{\circ }$ , the lift curves are no longer collapsed in the current measurements. At lower levels of turbulence, the wing exhibited traditional stall behaviour characterised by a sharp drop in lift between $\alpha = 17^{\circ }$ and $18^{\circ }$ . On the drag coefficient curve, $C_D = D / ((1/ 2) \rho U_{\infty }^{2} A )$ (where $D$ is the drag), for the low turbulent case shown in figure 5(b), this resulted in a sudden increase in drag. As the FST intensity was increased, the angle at which stall occurred shifted to larger $\alpha$ . This finding is consistent with several previous studies (Swalwell, Sheridan & Melbourne Reference Swalwell, Sheridan and Melbourne2001; Ahmadi-Baloutaki et al. Reference Ahmadi-Baloutaki, Carriveau and Ting2014; Kay et al. Reference Kay, Richards and Sharma2020; Thompson et al. Reference Thompson, Biler, Symon and Ganapathisubramani2023). The delay in stall with the addition of FST is likely a result of the boundary layer remaining attached at higher $\alpha$ and the suppression of the recirculation region, as demonstrated for similar FST conditions by Thompson et al. (Reference Thompson, Biler, Symon and Ganapathisubramani2023, Reference Thompson, Biler, Symon and Ganapathisubramani2025). Furthermore, no sudden stall was observed. Instead, the wing gradually lost lift while the drag increased up to $\alpha = 26^\circ$ .

Finally, to account for the fluctuations caused by the turbulence, the standard deviation of the lift and drag time series are shown in figures 5(c) and 5(d). These fluctuations increase with FST intensity in agreement with the findings of Thompson et al. (Reference Thompson, Biler, Symon and Ganapathisubramani2023). The observed increase in fluctuations suggests that FST has an impact on the instantaneous lift coefficient. Additionally, the drag curve displays a progressively increasing sensitivity to turbulence as $\alpha$ increases, with minimal fluctuations observed when approaching $\alpha = 0^{\circ }$ .

3.2. Wing deflections

In this section, the analysis of the wingtip displacements caused by the different turbulence cases measured with DIC (introduced in § 2.5) is presented. The statistical parameters, namely the mean and standard deviation of the deflections normalised with the wing’s maximum thickness, $\tau$ , characterising the wing’s deflection along the lift direction ( $\Delta y_w$ ) are presented in table 3. The maximum mean deflection of the wing is $12\,\%{-}13\,\%$ while the maximum standard deviation of this parameter is ${\approx}7\,\%$ . When normalised with the wing’s span instead of thickness, the mean and standard deviation of the deflections are ${\lesssim}1\,\%$ . Thus, the current wing model is certainly flexible, but far from the traditional ascription of the term aeroelastic. The mean of the wing’s (absolute) deflections and the corresponding standard deviations are seen to increase with FST. This indicates that the deflections are both greater and more violent as one increases the turbulence intensity.

Table 3. Mean of the wing’s (absolute) deflections and its standard deviation along the lift direction ( $\Delta y_w$ ) at the leading and trailing edges non-dimensionalised with the maximum thickness of the wing, $\tau$ .

Singular value decomposition (SVD) was applied on the fluctuating component of the DIC measurements to study the wing’s deformation in more detail, the results pertaining to which are presented in figure 6. The first SVD mode contains $13.2\,\%$ , $26\,\%$ , $92.5\,\%$ and $93.5\,\%$ of the total energy for cases REF, A, B and C, respectively, as shown in figure 6(a). These are the most energetic modes of the respective cases, with the second modes containing only ${\approx}1\,\%$ of the total energy.

The spatial organisation of the fluctuations corresponding to the most energetic mode for each tested case are presented in figure 6(b–e). Here, a coordinate system specific to the wing model $(x_w,y_w,z_w)$ is used. While $z_w$ is aligned with the $z$ direction of the main coordinate system, $x_w$ and $y_w$ are marginally different from the $x$ and $y$ coordinate owing to the $\alpha =5^{\circ }$ angle of the wing. The colour scale used here reflects relative values rather than absolutes, as the focus here is on the relative displacements of the wingtip points with respect to each other. These most energetic modes identify that the wingtip primarily bends along the lift direction (i.e. $\pm y$ ). Thus, the wing exhibits the behaviour of a cantilever beam. For cases REF and A, this mode is not particularly strong, while for cases B and C it is quite clear that the airfoil bends off at the tip. While it is not shown here, the second SVD modes corresponded to the twisting motion of the wing along its midchord. However, as previously stated, this mode contained an order-of-magnitude less energy than the first mode. This response of the wing’s structure is similar to the flexible wing experiments of Chen et al. (Reference Chen, Wang and Gursul2018b ) and Thompson et al. (Reference Thompson, Biler, Symon and Ganapathisubramani2025), where the bending deformation of the wingtip was shown to be much stronger than the torsional deformation with the addition of incoming fluctuations. Similar to the observations from the statistical analysis in table 3, SVD computations also indicate that the wing’s vibration increases with incoming FST. However, it seems to have a significant effect only on the bending motion of the wing and does not influence the wing’s twist along its midchord.

Figure 6. SVD of DIC measurements. The relative energy of the various modes compared with the total energy for the different cases is presented in (a). Spatial eigenfunction of the most energetic SVD mode of the different cases are shown in (b–d).

The frequency spectra of the displacement of the wingtip trailing edge are presented in figure 7(a). The natural frequency of the wing’s structure was measured in a static (no flow) environment by striking the airfoil and measuring its response with a laser displacement sensor (ILD2310-05) and was found to be $6.2\,\mathrm{Hz}$ . In all turbulence cases, a distinct energy peak is observed at the same frequency, which corresponds to a chord-based reduced frequency of $f c / U_{\infty } = 0.2$ . It is worth noting that while it is not presented here, the frequency content of the first SVD mode also corresponds to the same value. Thus, the magnitude of the wing’s vibration is sensitive to the level of incoming turbulence, but the specific active frequency still remains the natural frequency of the wing.

Finally, as mentioned in § 3.1, the lift force measured using the force balance showed an increase in lift fluctuations with incoming turbulence intensity. This is especially apparent for cases B and C in figure 5(c). The spectra of these lift fluctuations when the wing is at $\alpha = 5^{\circ }$ are presented in figure 7(b). The energy in the spectra of the cases REF and A are very low, which corroborates the observations in figure 5(a). Interestingly, however, for the high turbulence cases B and C, the dominant frequency in the lift fluctuations is observed to also be at $6.2\,\mathrm{Hz}$ . Thus, when the amplitude of the wing’s vibration is high enough, the generated lift force is not constant any more, but fluctuates at the same frequency as the wing’s vibrations.

Figure 7. (a) Power spectra of the wingtip trailing edge displacement. The magnitude of the free vibration test spectrum was shifted vertically to ease comparison with the other spectra. (b) Power spectral density (PSD) computed on lift fluctuations measured using the force balance for $\alpha = 5{^{\circ }}$ . Spectra of cases A, B and C are shifted by 1, 2 and 3 decades, respectively, for clarity.

4.1. Swirling strength and vortex size

The swirling strength ( $\lambda _{ci}$ ), which is the imaginary part of the complex eigenvalue of the local velocity gradient tensor, is used to identify and depict the wingtip vortex. The contour maps of the time-averaged streamwise swirling strength, $\overline {\lambda }_{ci}$ , for each test case are presented in figure 8(a–d). Here, the maximum swirling strength at the centre of the vortex was used to non-dimensionalise the rest of the field in order to maintain a consistent colour scale for all cases. The maximum values of the time-averaged swirling strength, normalised with global parameters, i.e. $\mathrm{max}(\overline {\lambda }_{ci}) \boldsymbol{\cdot }c^2/U_{\infty }^2$ , decreased with FST as shown in figure 8(i). In figure 8(a–d), a high concentration of swirling strength that is circular in shape is observed near the centre of the frame. The size of the swirling strength contours appear to grow with increasing FST. As highlighted in § 1, previous studies have demonstrated that an increase in the incoming turbulence intensity leads to an increase in the meandering amplitude of the wingtip vortex (Bailey & Tavoularis Reference Bailey and Tavoularis2008; Pentelow Reference Pentelow2014; Ben Miloud et al. Reference Ben Miloud, Dghim, Fellouah and Ferchichi2020). In this context, the time-averaged swirling strength in figure 8(a–d) would naturally give the impression of a vortex that occupies more area. In reality, this figure does not provide information about the diffusion of the vortex, but is an indicator of how the FST increases vortex meandering.

Figure 8. (a–d) Normalised global time-average of the swirling strength ( $\overline {\lambda }_{ci}$ ) of the four test cases (no alignment of vortex-centres). The black dashed line represents the surface projection of the wing. (e–h) Normalised, conditionally time-averaged swirling strength ( $\hat {\lambda }_{ci}$ ) of the four test cases (alignment of vortex centres). (i) Maxima of the time-averaged and conditionally time-averaged swirling strengths plotted against $u'/U_{\infty }$ . Note that the axes of (e–h) are smaller than (a–d).

To visualise and study the changes caused by the FST on the wingtip vortex itself, the swirling strength data was averaged across time after adjusting each instantaneous field to be such that the origin is at the vortex centre (i.e. the vortex-based coordinate system $(\hat {y}, \hat {z}) = (0,0)$ ), following the data reduction method of Heyes et al. (Reference Heyes, Jones and Smith2004). This is henceforth referred to as conditional (time-)averaging (also introduced at end of § 2.1). To identify the centre of the vortex in each frame, a non-Galilean invariant approach introduced by Graftieaux, Michard & Grosjean (Reference Graftieaux, Michard and Grosjean2001) was implemented. This approach consists of computing the $\varGamma _{1}$ criterion from the SPIV fields to extract the position of the vortex centre, with $\varGamma _{1} \in [0,1]$ defined as a scalar function,

(4.1) \begin{equation} \varGamma _{1}(P) = \frac {1}{N} \sum _{S} \frac {\left (\boldsymbol{R}_{\textit{PM}} \times \boldsymbol{U}_{M}\right ) \boldsymbol{\cdot }\boldsymbol{e}_{x}}{\left |\boldsymbol{R}_{\textit{PM}}\right |\left |\boldsymbol{U}_{M}\right |} = \frac {1}{N} \sum _{S} \operatorname {sin}\left (\theta _{M}\right)\!, \end{equation}

where $S$ is a 2-D area surrounding a fixed point defined by $P$ , $M$ is any point lying in $S$ , $\boldsymbol{e}_{x}$ is the unit vector normal to the measurement plane and $\theta _{M}$ represents the angle between the velocity vector, $\boldsymbol{U}_{M}$ , and the radius vector, $\boldsymbol{R}_{\textit{PM}}$ . When the angles $\theta _{M}$ between the velocity vectors at points $M$ and the radius vector approach $\unicode{x03C0} /2$ for a point $P$ , $\varGamma _{1}$ tends towards unity, the maximum value for $\varGamma _{1}$ . In practice, $\varGamma _{1}$ rarely reaches unity and the point of maximum $\varGamma _{1}$ is considered the vortex centre.

The result of this procedure (i.e. determining the instantaneous vortex positions and recentring each swirl strength field before averaging across frames) is shown in figure 8(e–h). In these plots, the conditionally averaged swirl strength ( $\hat {\lambda }_{ci}$ ) is normalised by the maximum value in the respective field. This normalisation highlights the relative distribution of swirling strength around the vortex core, which appears to broaden slightly with increasing FST, suggesting a slight diffusion of the vortex. Nonetheless, in this frame the change to the vortex is significantly smaller than the laboratory-frame averaging suggests. The corresponding maximum values of the conditionally averaged swirl strength for the different cases normalised with global parameters are presented in figure 8(i). A linear decrease of this parameter is observed with increasing FST, indicating a reduction in vortex strength, consistent with previous findings (Bailey et al. Reference Bailey, Tavoularis and Lee2006; Ghimire & Bailey Reference Ghimire and Bailey2017; Ben Miloud et al. Reference Ben Miloud, Dghim, Fellouah and Ferchichi2020).

Figure 9. (a-d) Evolution of vortex core radius versus normalised time for the four test cases. Dashed (white) line represents mean radius. (e) Probability density function (PDF) of the vortex size.

The results in figure 8(e–h) indicate that the average vortex size remains fairly constant with FST. To investigate this further, figure 9(a–d) display the time evolution of the vortex-core radius, with the corresponding average value $r_{c} (=\langle \hat {r}_{c} \rangle )$ marked by dashed white lines. Herein, the vortex core radius is defined as the radial position of the maximum azimuthal velocity (see Appendix A). The vortex size remained relatively stable for cases REF and A. However, for case A, the deviation from the mean value increased, suggesting a slight decrease in the overall stability. For the more turbulent cases B and C, the deviation from the mean increases further, indicating larger changes to the size of the vortex in time. The average vortex radius, $r_{c}$ , is shown through the PDF in figure 9(d), where a minor increase in average vortex size can be observed, with $r_{c}/c=0.0743$ for case REF increasing to $r_{c}/c=0.0828$ for case C. This is paired with an increase in scatter and a slight asymmetry in favour of larger vortices.

It is worth noting that azimuthal velocity profiles and the circulation within the tip vortex were also studied. These analyses are not included here for conciseness but can be found in Appendix A and are in good agreement with previous studies (Bailey et al. Reference Bailey, Tavoularis and Lee2006; Pentelow Reference Pentelow2014; Ghimire & Bailey Reference Ghimire and Bailey2017; Ben Miloud et al. Reference Ben Miloud, Dghim, Fellouah and Ferchichi2020). Notably, the maximum azimuthal velocity within the vortex core decreases with increasing FST and it is also inferred that FST enhanced the rate at which the vortex reached its fully developed structure (Bailey et al. Reference Bailey, Tavoularis and Lee2006; Ghimire & Bailey Reference Ghimire and Bailey2017). A possible reduction in circulation near the core radius indicating an increase in vortex diffusion with FST was also observed, which is in agreement with Ghimire & Bailey (Reference Ghimire and Bailey2017) and Ben Miloud et al. (Reference Ben Miloud, Dghim, Fellouah and Ferchichi2020). This observation, in conjunction with a decrease in the maximum streamwise swirl strength (figure 8 i) and a slight increase in the average radius of the tip vortex (figure 9), support the observation made with figure 8(e–h) that the diffusion of the vortex increases with FST.

4.2. Probability distributions of vortex position and meandering amplitudes

Bailey & Tavoularis (Reference Bailey and Tavoularis2008) and Van Jaarsveld et al. (Reference Van Jaarsveld, Holten, Elsenaar, Trieling and Van Heijst2011) have previously reported that the meandering amplitude of the vortex, denoted in the current work by $\sigma _i$ ( $i$ representing coordinate direction), increased with FST. Figure 10 presents the 2-D probability density distribution of the instantaneous vortex centres relative to its mean position. It is worth noting that it is this conditioned (fluctuations) signal that is used in the rest of this section when referring to vortex meandering. The darker colours in figure 10 represent the areas where the vortex centre was more likely to be located. For REF, the results show a highly localised vortex centre that remained stationary in time. For case A, the circular concentration area of the probability of the vortex centre location starts to extend, although it still remained fairly localised. In the two most turbulent cases (B and C), the PDFs indicate that the instantaneous vortex centre could be located in a wider area. Thus, as FST increased, there is a noticeable increase in the scattering of the vortex positions, in agreement with Bailey & Tavoularis (Reference Bailey and Tavoularis2008) and Van Jaarsveld et al. (Reference Van Jaarsveld, Holten, Elsenaar, Trieling and Van Heijst2011), indicating an increase in vortex meandering.

Figure 10. Probability density function of the instantaneous vortex centre relative to the mean vortex centre position (meandering conditioned) of the four different cases. A least-squares-based ellipse fit with its minor and major axis (dashed and dash–dotted, respectively) is also shown.

To determine the meandering amplitudes in the transverse and vertical directions, denoted by $\sigma _{y}$ and $\sigma _{z}$ , respectively, the standard deviations of the time series for the vortex axis displacements from the mean position were calculated for each direction. The dependence of these amplitudes on the incoming turbulence intensity is presented in figure 11. As expected and noted by Bailey et al. (Reference Bailey, Pentelow, Ghimire, Estejab, Green and Tavoularis2018) and Ben Miloud et al. (Reference Ben Miloud, Dghim, Fellouah and Ferchichi2020), the meandering amplitudes strongly increase in both directions with FST. While previous efforts demonstrated this only for turbulence intensities less than $6\,\%$ , the current results show that this tendency remains true even at the highest tested turbulence intensity of $13\,\%$ . Furthermore, the meandering is homogeneous in both directions of the transverse plane.

Figure 11. Variation of meandering amplitude ( $\sigma _i$ ) with incoming turbulence intensity: $*$ , present results; $\triangle$ , Bailey & Tavoularis (Reference Bailey and Tavoularis2008) ( $x/c = 3$ , ${\textit{Re}}_{c} = 2.4 \times 10^{5}$ ); $\lozenge$ , Pentelow (Reference Pentelow2014) ( $x/c = 2.85$ , ${\textit{Re}}_{c} = 2.4 \times 10^{4}$ ); $\circ$ , Ben Miloud et al. (Reference Ben Miloud, Dghim, Fellouah and Ferchichi2020) ( $x/c = 2.5$ , ${\textit{Re}}_{c} = 2.0 \times 10^{5}$ ); $\triangleleft$ , Bölle et al. (Reference Bölle, Brion, Couliou and Molton2023) ( $x/c = 2$ , ${\textit{Re}}_{c} = 1.7 \times 10^{5}$ ).

Finally, ellipses were fitted to the vortex-position data using a least-squares method, and are plotted on the PDFs in figure 10. For this fit, only the data points where the tip vortex was within the region defined by $\hat {y}/c = \pm \sigma _y \boldsymbol{\cdot }\sqrt {2}$ and $\hat {z}/c = \pm \sigma _z \boldsymbol{\cdot }\sqrt {2}$ were considered. The ratio of the major and minor axes of the ellipse was ${\approx} 1.1$ for all cases, which suggests that the meandering motion was approximately homogeneous in the transverse plane. The orientation or tilt of the ellipse, however, increased with FST. It was negligible for case A but increased to $0.06\unicode{x03C0}$ and $0.2\unicode{x03C0}$ for cases B and C, respectively. This is an indication that the meandering tends in the direction of the lift force with increasing FST intensity. In the flexible wing experiments of Chen et al. (Reference Chen, Wang and Gursul2018b ), probability distribution of the instantaneous vortex location was more dispersed downstream of the flexible wing compared with their rigid wing. This suggests that the current observations may be attributed to the wing’s deflections, a point that will be explored further in § 5.

4.3. Spectral signature of wingtip-vortex’ motion

The frequency content of the vortex core’s motion is investigated further, to complement the results of the meandering amplitude. Given that the dominant vibration of the wing is along the lift (i.e. $y$ ) direction as shown in § 3.2, only the meandering spectra associated with the transverse direction, $\phi _{yy}(f)$ , are presented and discussed further.

Figure 12. Premultiplied frequency spectra of the vortex motion $\phi _{yy}$ . Dash–dotted (red) line indicates normalised frequency of the wing’s vibration. Note that the vertical axes of (a) and (b) are two orders-of-magnitude lower than (c) and (d).

The premultiplied frequency spectra of the vortex motion for each turbulent case, as presented in the literature (see Bailey et al. Reference Bailey, Pentelow, Ghimire, Estejab, Green and Tavoularis2018; Chen et al. Reference Chen, Wang and Gursul2018a ) are shown in figure 12. The vertical (red) dash–dotted lines correspond to the wing’s vibration frequency ( $6.2\,\mathrm{Hz}$ ) in normalised form. Note that, in some frames, the identification of the vortex centre using the $\varGamma _{1}$ criterion was not conclusive; often because, for the higher turbulence conditions, the vortex moved out of SPIV FOV for short periods. Consequently, the amount of information lost in these cases was $0\,\%$ , $0\,\%$ , ${\approx}1\,\%$ and ${\approx}10\,\%$ for cases REF, A, B and C, respectively. The fillgaps () function in MATLAB, which estimates the missing values through forward and backward autoregressive fits of the remaining data, was employed to fill the missing data. This was validated by taking complete sections of the data and removing data points and refilling them with this approach; the statistics and spectra were found to be preserved.

Figure 10 demonstrates that vortex meandering motion for the REF case is minimal. Thus, the spectrum in figure 12(a) shows no specific dominant frequency. With the addition of FST (figure 12 b–d), the meandering energy is seen to be in a frequency band spread over two orders-of-magnitude. Between the least (A) and most turbulent (C) cases, the energy under the curve has a difference of two orders-of-magnitude. These observations are consistent with the results shown earlier regarding vortex position in the $y{-}z$ plane (figure 10) and the amplitude of meandering (figure 11).

The broadband nature of the meandering spectra centred around some frequency band was shown also by Bailey et al. (Reference Bailey, Pentelow, Ghimire, Estejab, Green and Tavoularis2018). This demonstrates the stochastic nature of tip-vortex meandering and that the measured meandering is composed of displacements at different scales. The fluctuations in this broad frequency band correspond to relatively long wavelengths ( $\lambda = U_{\infty }/f$ ) in the range $0.74\lt \lambda /c \lt 83$ ( $10 \lt \lambda /r_c \lt 1000$ ). It has been shown previously that even in low turbulence conditions ( $u'/U_{\infty }\approx 0.5\,\%$ ), the long wavelength motions are more dominant than short wavelengths (Chen et al. Reference Chen, Wang and Gursul2018a ). Furthermore, based on the observations of Bailey et al. (Reference Bailey, Pentelow, Ghimire, Estejab, Green and Tavoularis2018), the addition of FST mostly dampens short-wavelength motions. Considering both these observations in the previous literature, the short-wavelength vortex motion is not expected to be dominant in the current results either.

Together with the energy within the frequency band associated with vortex motion increasing, a reduction in the peak frequency (both normalised and dimensional) is also seen in figure 12(b–d). One can thus identify that the strongest vortex motion occurs at lower frequencies (or larger meandering wavelengths) as the intensity of incoming turbulence increases. A similar observation was made by Bailey et al. (Reference Bailey, Pentelow, Ghimire, Estejab, Green and Tavoularis2018) when comparing their low-turbulence case and high-turbulence cases in the near wake. However, the current set-up includes a flexible wing which could have additional effects on the vortex dynamics on top of the incoming turbulence.

Artefacts of the wing’s vibrations are noticeable in the spectra of all cases with the active grid. In case A (figure 12 b), a small peak is noticed at the normalised frequency of $0.2$ which could be attributed to the vibration of the wing, but it is an order-of-magnitude lower in energy compared with the peak in the spectrum. A similar observation can be made in figure 12(c) showing the spectrum of case B. The spectral peak at the normalised frequency of the wing’s vibration is now close to the maxima of the frequency band. In the most turbulent case C in figure 12(d), the wing’s vibration frequency is within the dominant vortex motion’s frequency range. Thus, these results strongly suggest that the spectral content of the vortex core’s displacement is linked to the frequency of the wingtip motion.

4.4. Proper orthogonal decomposition assessment of vortex meandering

More insight into the dynamics of wingtip vortices can be obtained using POD (see, for example, Edstrand et al. Reference Edstrand, Davis, Schmid, Taira and Cattafesta III2016, Chen et al. Reference Chen, Wang and Gursul2018a , Reference Chen, Wang and Gursulb and Dghim et al. Reference Dghim, Ben Miloud, Ferchichi and Fellouah2021). The same is done here to further investigate the impact of FST on the wingtip vortex. The snapshot formulation of POD proposed by Sirovich (Reference Sirovich1987) is employed for this purpose. It was applied only on the fluctuation fields of the transverse velocity components. In doing so, the fluctuating velocity fields arranged in matrix form as $Q=[(v-\overline {v}) ; (w-\overline {w})]$ are decomposed in the following manner:

(4.2) \begin{equation} { Q({y},{z},t) = \sum _{n} \varLambda _n \varTheta _n(t) \varPsi _n({y},{z}). } \end{equation}

Here, $\varPsi _n({y},{z})$ are a set of empirical eigenfunctions (spatial POD modes) that contain information regarding the spatial organisation of various fluctuations with corresponding eigenvalues, $\varLambda _n$ , representing the modal energy and temporal coefficients, $\varTheta _n(t)$ , with $n$ being the mode number. It is worth noting that the modes are sorted in descending order of their relative modal energy. Henceforth, any mention of mode number corresponds to its rank based on its relative modal energy. Furthermore, given that the amplitude of meandering was very low for the REF case (figure 11) and that the frequency content of the vortex’ motion in this flow condition is not well resolved in the current data as seen in figure 12(a), it is not considered in the current analysis. The first four spatial modes ( $y$ component) of cases A, B and C are shown in figure 13 as filled contours. The cross product of the two components of the POD spatial modes, which can be related to vorticity along the streamwise direction ( $\varOmega _n^x = \boldsymbol{\nabla }\times \varPsi _n(y,z)$ ), is computed and also shown in these plots as green and pink contour lines. The relative modal energy of the first 10 eigenmodes is presented in figure 14(a). The spectra of the temporal coefficients of the first four POD modes are plotted in figure 14(b–d) for cases A, B and C, respectively.

Figure 13. Spatial organisation of fluctuating fields corresponding to POD eigenmodes $n = 1$ (a–c), $n=2$ (d–f), $n=3$ (g–i) and $n=4$ (j–l) for cases A (a,d,g,j), B (b,e,h,k) and C (c,f,i,l). Filled contours represent the component of POD spatial modes along the lift direction, $\varPsi _n^{{y}}$ (blue, positive; red, negative). Contour lines represent the vorticity-like term computed with both components of the respective POD spatial modes ( $\varOmega _n^x = \boldsymbol{\nabla }\times \varPsi _n({y},{z})$ ; green, positive; pink, negative).

When looking at the $\varOmega _n^x$ contour lines in figure 13, the first modes of all three cases contain two regions of oppositely signed vorticity, forming a vortex dipole together. A similar organisation of fluctuations is observed in the second mode of cases A and B as well. In case C (figure 13 f), while there are two structures observed, the shape of the one on the right-hand side is not as well-resolved as the rest. Nevertheless, these appear to be paired in all three cases and orthogonally rotate together. As such, these mode-pairs could be associated with helical displacement modes arising from a Kelvin wave with an azimuthal wavenumber of $m = 1$ as was conjectured by previous studies as well (Fabre, Sipp & Jacquin Reference Fabre, Sipp and Jacquin2006; Dghim et al. Reference Dghim, Ben Miloud, Ferchichi and Fellouah2021). A deeper analysis of the vortex stability would be useful to confirm this in the future. Nevertheless, this dipole is a result of the transverse velocity components being non-zero at its centre, which leads to vortex meandering (Fabre et al. Reference Fabre, Sipp and Jacquin2006; Edstrand et al. Reference Edstrand, Davis, Schmid, Taira and Cattafesta III2016; Chen et al. Reference Chen, Wang and Gursul2018a ,Reference Chen, Wang and Gursul b ; Dghim et al. Reference Dghim, Ben Miloud, Ferchichi and Fellouah2021). This last idea is in fact supported by the POD results presented here. The filled contours in figure 13(a–f), which represent the organisation of transverse velocity fluctuations in the $y$ direction ( $\varPsi _n^y(y,z)$ ), is non-zero at the centre of the dipole in the first two modes of all three cases.

The structures in these spatial modes also increase in size with FST, which is synonymous with the increase in the size of the swirl-strength contour in figures 8(a–d). The relative energies of these mode-pairs also increase with the addition of turbulence as seen in figure 14(a). For case A, the first mode-pair contains ${\approx} 17\,\%$ of the total energy, which increases to approximately $41\,\%$ for cases B and C. Given that the spatial modes of this mode-pair resemble helical displacement modes previously linked to wingtip vortex meandering (Edstrand et al. Reference Edstrand, Davis, Schmid, Taira and Cattafesta III2016; Chen et al. Reference Chen, Wang and Gursul2018a ,Reference Chen, Wang and Gursul b ; Dghim et al. Reference Dghim, Ben Miloud, Ferchichi and Fellouah2021), these two features observed in the POD analysis point towards the increase in meandering amplitude with the addition of FST discussed in § 4.2. The premultiplied spectra of the temporal coefficients of these modes in figure 14(b–d) contain a spectral signature comparable to that of the meandering motion shown in figure 12 as well.

Figure 14. (a) Relative energies of the first 10 POD eigenmodes of all cases. (b–d) Premultiplied spectra of temporal coefficients of first four POD modes for cases A, B and C, respectively. Spectra of various modes are weighted with their relative energy for clarity. Dash–dotted line indicates normalised frequency of wing’s vibration.

The third and fourth spatial modes of case A (figures 13 g and 13 j) do not show any clear structures. Their relative energies have also dropped below $5\,\%$ . On the contrary, the third and fourth most energetic modes of cases B and C are vortex deformation modes. Specifically, the spatial organisation of the fourth mode resembles those resulting from the amplification of elliptical instabilities within the vortex tube with an azimuthal wavenumber of $m=2$ (see Lacaze, Ryan & Le Dizes Reference Lacaze, Ryan and Le Dizes2007). These modes represent only approximately $12\,\%$ of the total energy which is less than a third of the first two displacement modes. While the wing’s increased vibration and/or the turbulence might be amplifying these deformation modes, the meandering motion is far stronger.

Interestingly, the orientation of structures in the most energetic modes across all three cases (figure 13 a–c) suggests vortex displacement along the $y$ direction, which coincides with the direction of the wingtip’s deflection when FST is present. While this mode is only slightly more energetic than the second ( ${\approx}5\,\%$ ), this subtle preference in vortex displacement may suggest a potential influence of the wing’s vibrations on vortex motion.

5.1. Wing and vortex motion

The relationship between the wing’s deflection and vortex meandering is investigated here. To do so adequately, a common reference signal, for both the motion of the wing and the tip vortex, is necessary. However, such a signal is not available in the current experimental set-up as both the wing and vortex were allowed to freely move with the incoming turbulence. The Hilbert transform, which is the convolution of a real valued signal with the Cauchy kernel ( $1/\unicode{x03C0} t$ ), helps circumvent this issue with the Cauchy kernel playing the role of the common reference. In fact, this method was recently employed by Kushwaha et al. (Reference Kushwaha, Worth, Dawson, Gupta and Li2022) to assess the phase difference between the shear layers rolling up on either side of the potential core of a round jet subjected to transverse acoustic forcing. In the current study, the Hilbert transform was applied to the time series of the wingtip trailing edge displacement obtained from DIC and that of the conditioned vortex meandering from SPIV, which provided the instantaneous phase of the two motions along the $y$ direction ( $\theta _{yy,W}$ and $\theta _{yy,V}$ , respectively) with respect to Cauchy kernel. The phase difference between the two signals was then computed and plotted as a probability distribution presented in figure 15.

Figure 15. Probability density function of the instantaneous phase difference between the wing’s deflection and vortex meandering motion. Here $\theta _{yy,W}$ and $\theta _{yy,V}$ are the instantaneous phases of the wing and vortex motion, respectively, obtained through the Hilbert transform.

In cases REF and A, the probability distributions show no dominant phase relation, with the occurrence of all phase differences having a probability of ${\approx}0.5$ . For cases B and C, a dominant phase difference centred around $0.5\unicode{x03C0}$ is evident, indicating a relation between the wing and vortex motions. Specifically, it shows that the wing’s deflection leads that of the vortex meandering motion by $0.5\unicode{x03C0}$ . This finding, that a dominant phase relation between the wing and vortex motions appears with FST, can be attributed to two already noted observations. First, the vortex motion moves towards longer wavelengths or smaller frequencies with FST, as discussed in §§ 4.3 and 4.4, which in the current experimental set-up also means a drift towards the wing’s vibration frequency. Second, the wing’s vibration amplitude increases significantly between cases REF and A, and cases B and C, as seen in figure 7(a). The results of this instantaneous phase-relation analysis follows the same trend, suggesting that the meandering motion of the vortex axis is correlated with the displacement of the wing at higher FST.

5.2. Effect on the wingtip vortex

It is imperative to also investigate the effect of the wing’s deflection on the wingtip vortex itself in order to fully characterise its effect on the wake. The Hilbert transform analysis in § 5.1 shows that the wing’s deflections has an effect on the vortex’ motion. A similar trend was observed in the POD analysis of the full SPIV FOV presented in § 4.4. In that analysis, a substantial portion of the total energy was concentrated in the first mode pair, representing vortex meandering ( ${\approx} 17\,\%$ for case A, ${\approx} 41\,\%$ for cases B and C). In contrast, the third and fourth most energetic modes in cases B and C, resembling deformation modes, accounted for only ${\approx} 12\,\%$ of the total energy. To better isolate the effects on the vortex itself, it is necessary to remove the vortex core’s meandering motion before performing POD, which is undertaken here. To this effect, velocity vectors within a square window of size $\pm 2r_c$ centred around the vortex centre were extracted from the original full-field time series. It is worth noting that the meandering component was not fully removed and the vortex did have some motion within the $\pm 2r_c$ interrogation window. Nevertheless, snapshot-POD was applied on this extracted time series of cases A, B and C, and only on the fluctuations of the transverse velocity components. Similar to the previous POD computation, the fluctuating velocity fields arranged in matrix form as ( ${Q}=[(v-\hat {v}) ; (w-\hat {w})]$ ) are decomposed using (4.2). Once again, $\varPsi _n(\hat {y},\hat {z})$ are a set of empirical eigenfunctions (spatial POD modes) that contain information regarding the spatial organisation of various fluctuations with corresponding eigenvalues, $\varLambda _n$ , representing the modal energy and temporal coefficients, $\varTheta _n(t)$ , with $n$ being the mode number. It is worth noting that the modes are sorted in descending order of their relative modal energy here as well.

Figure 16. (a) Relative energies of the first $10$ POD eigenmodes of all cases. (b–d) Premultiplied spectra of temporal coefficients of first three POD modes for cases A, B and C, respectively. Spectra of various modes are weighted with their relative energy for clarity. Dash–dotted line indicates normalised frequency of wing’s vibration.

The relative energy of the first 10 most energetic POD modes corresponding to the three cases are presented in figure 16(a). In case A, the energy of subsequent modes decreases monotonically, with the first mode containing ${\approx} 10\,\%$ of the total energy. For test cases B and C, the energy of the first three POD modes is much higher compared with those in case A. The first modes contain ${\approx} 25\,\%$ and ${\approx} 30\,\%$ of the total energy in cases B and C, respectively, which is significantly higher compared with the first mode of case A. The results reported here (see § 4) and previous studies (Ghimire & Bailey Reference Ghimire and Bailey2017; Bailey et al. Reference Bailey, Pentelow, Ghimire, Estejab, Green and Tavoularis2018; Ben Miloud et al. Reference Ben Miloud, Dghim, Fellouah and Ferchichi2020) have demonstrated that increasing FST leads to a reduction in the vortex strength. Thus, the increase in the relative energy of the first three POD modes of cases B and C, while that of the other modes remaining almost constant (less than $1\,\%$ change from case A), cannot be due to the FST’s action on the tip vortex. It can also not be due to the meandering motion as this has effectively been removed through the windowing around the vortex centre. This is rather attributed to the increased amplitude of the wing’s deflection as will be discussed later.

Premultiplied spectra of the temporal coefficients of the first three POD modes corresponding to the three cases are shown in figure 16(b–d). Similar to figure 14(b–d), the spectra of individual modes are weighted by their relative energy to improve clarity. No dominant peaks are observed in the spectrum of case A (figure 16 b). With further increase in FST, the frequency band with the maximum energy tends towards the vibration frequency of the wing, with the frequency peak in the spectrum almost at the wing’s vibration frequency for case C (figure 16 d). This spectral signature is expected to be from the deflections of the wing, indicating the effect it has on the vortex itself, separate from how the vortex meanders.

Figure 17. Spatial organisation of fluctuating fields corresponding to POD eigenmodes $n = 1$ (a–c) and $n=2$ (d–f) for cases A (a,d), B (b,e) and C (c,f). Filled contours represent component of POD spatial modes along the lift direction, $\varPsi _n^{\hat {y}}$ (blue, positive; red, negative). Contour lines represent the vorticity-like term computed with both components of the respective POD spatial modes ( $\varOmega _n^x = \boldsymbol{\nabla }\times \varPsi _n(\hat {y},\hat {z})$ ; green, positive; pink, negative).

The spatial organisation of fluctuations are now looked into to elucidate these effects further. The $\hat {y}$ component of the spatial eigenfunctions of the first two POD modes for cases A, B and C are shown in figure 17 as filled contours. The modes of case A do not contain a distinct structure. This is in agreement with the spectra of the corresponding temporal coefficients not indicating a dominant frequency, but also with the full-field POD computation of this case where no clear structures were visible in the third most energetic mode (figure 13 g) . In comparison, the first mode of cases B and C have a distinct pattern where the bottom-half is oriented along $+\hat {y}$ and the top is along $-\hat {y}$ . Clear structures are also noticeable in the second modes of these two cases.

The cross product of the two components of the POD spatial modes ( $\varOmega _n^x = \boldsymbol{\nabla }\times \varPsi _n(\hat {y},\hat {z})$ ), is computed and also shown in figure 17 as green and pink contour lines. In doing so, a distinct structure arises which indicates that the fluctuations represented by the first modes of cases B and C are organised in such a manner that it results in a rotational motion. This structure is similar to the elliptical structure in the third mode of cases B and C in the previous POD computation (figures 13 h and 13 i). In the second POD modes, two smaller rotational structures become evident through this process.

The motion of the vortex in the transverse plane or vortex meandering has been negated in the POD computation as previously explained. Thus, the fluctuations observed in the temporal coefficients’ spectra in figure 16 cannot represent any vortex motion in the transverse plane. Instead, they are associated with the strength (and size) of the rotational structures in figure 17 (especially mode 1 of cases B and C) varying in time. Evidence of this can be found in figure 18, where the normalised streamwise vorticity of six instantaneous velocity fields spread along one cycle of the wing’s vibration of $6.2\,\mathrm{Hz}$ is presented for the most turbulent case C. These instantaneous fields were reconstructed with the most energetic POD eigenmodes that account for $95\,\%$ of the total energy. The clear variation of the vorticity level and small changes to the vortex size with time are identifiable in these snapshots.

To get a global picture of this phenomenon, the instantaneous vortex circulation was computed from the reconstructed fields of all three turbulence cases within a square window of size $\pm r_c$ around the vortex centre. The premultiplied spectra of this quantity are presented in figure 19. The spectral content in case A is extremely low as expected, given the similar observations in figure 16(b). For cases B and C, fluctuations with maximum spectral content tend towards the wing’s normalised vibration frequency, with it being almost same for case C, as was observed with the spectra of the vortex core’s motion as well. This is further proof of the temporal variation of the wingtip vortex’ strength in the current experimental study. Having demonstrated that these fluctuations, if not originating from the wing’s deflection itself, are at least correlated with it, it can be safely said that the variation in the strength of the tip vortex in time must be a result of the wing’s vibration.

Figure 18. Normalised streamwise vorticity computed on six instantaneous velocity fields reconstructed from POD eigenmodes of case C. The six instantaneous fields represent one cycle of the wing’s vibration at $6.2\,\mathrm{Hz}$ .

Figure 19. Premultiplied spectra of the circulation $\varGamma (t)$ computed from the reconstructed instantaneous velocity fields for the three turbulent cases. Dash–dotted (red) line indicates normalised frequency of wing’s vibration. Note that the spectral content in case A is extremely low in comparison with cases B and C.

This observation can be explained as such. In the case of a rigid wing, the lift force generated by the wing and that perceived by the vortex is the same. Thus, a temporal variation in the vortex strength does not occur. This is in fact the observation with case A, where the vibration amplitude of the wing is at least three orders-of-magnitude lower than cases B and C as shown in figure 7(a) through the DIC measurements. The same does not hold when the vibration amplitude of the wing increases. In this scenario, when the wingtip is moving in the $+y$ direction or along the lift, the cumulative lift force that is perceived by the flow is in fact lower. When it is moving in the direction opposite to the lift (i.e. along $-y$ ), the cumulative force is higher. In fact, this temporal variation in the lift is represented in figure 5(a) by the vertical bars, which, as shown in figure 7(b), has a spectral peak at the wing’s vibration frequency as well. This fluctuation in the perceived lift force must lead to a variation in the strength of the wingtip vortex but also its radius as seen in figure 9.