Biomimetic Groove and Elliptical Bluffness Synergy for Enhanced Vortex-Induced Vibration Excitation

1  Introduction

The accelerating global transition toward renewable and distributed energy systems has prompted significant research interest in extracting useful power from fluid-structure interaction phenomena, particularly from low-speed hydrokinetic currents in rivers, tidal channels, and shallow coastal ocean streams. While large-scale hydraulic turbines remain the dominant technology for hydrokinetic conversion, there exists a rapidly expanding demand for compact, maintenance-free micro-energy harvesting systems designed to power autonomous Internet of Things (IoT) sensors, deep-ocean instrumentation, and remote wireless structural health monitoring nodes. In such deployment environments, including subsea pipelines, riverbed monitoring stations, and internal flow conduits, the available flow velocities are inherently extremely low (U∞<0.1 m/s), confining the operating Reynolds number strictly within the laminar regime (100≤Re≤500) [1,2]. Among the fluid-structure interaction mechanisms available at these ultra-low flow speeds, Vortex-Induced Vibration (VIV) has emerged as one of the most promising physical principles for small-scale energy conversion [3,4]. As a fluid flows past a bluff body, alternating low-pressure vortices detach periodically from its upper and lower surfaces, forming a characteristic pattern known as the Kármán vortex street. These alternating vortices generate a periodic cross-flow lift force oscillating at the Strouhal frequency, which drives an elastically mounted structure into sustained transverse oscillations [5,6]. When the vortex shedding frequency approaches the natural frequency of the mechanical system, a condition termed lock-in or synchronization, oscillation amplitudes are substantially amplified, providing the large-displacement, sustained mechanical motion required for efficient energy extraction [7,8]. The fundamental engineering appeal of VIV energy harvesting lies in its entirely passive, self-sustaining character: no active control systems, no rotating machinery, and no external power input are required. The seminal VIVACE (Vortex Induced Vibration Aquatic Clean Energy) concept introduced by Bernitsas et al. [1] demonstrated that properly designed elastically mounted cylinders can achieve lock-in oscillations in marine currents as slow as 0.25 m/s, establishing the practical viability of the VIV paradigm. Investigating these mechanisms at low Reynolds numbers therefore constitutes not merely a theoretical computational simplification, but a highly realistic operating condition for a well-defined and growing class of next-generation micro-power generation devices.

Nature has long provided engineering inspiration for passive aerodynamic control strategies. Among the most extensively studied biomimetic surface modifications are circumferential longitudinal grooves and surface dimples, whose geometric principles draw directly from biological archetypes such as the ribbed epidermis of the cactus (Cereus repandus) and the sinusoidal leading-edge tubercles of the humpback whale (Megaptera novaeangliae) flipper [9,10]. The cactus spine-groove pattern has proven particularly effective as a passive boundary layer manipulation strategy: the convex ridges separating adjacent grooves act as localized tripping elements that promote turbulent mixing in the near-wall shear layer, effectively delaying the onset of flow separation and reducing the extent of the near-wake recirculation region [11]. The direct aerodynamic consequence is a substantial reduction in the mean pressure drag, a narrowing of the wake, and, crucially from the perspective of structural vibration, a suppression of the alternating lift force amplitude [12,13]. The numerical and experimental investigations have consistently confirmed these beneficial stabilizing effects: Babu and Mahesh [14] performed direct numerical simulations of flow past cactus-shaped cylinders at Re=20, 100, and 300, reporting drag reductions of up to 22% and lift coefficient reductions exceeding 47%, attributed to quiescent recirculating flow within the cactus cavities. Furthermore, Talley et al. [15] reported analogous drag reduction behavior across a range of groove wave numbers. The overwhelming consensus in the published literature therefore positions biomimetic surface grooves exclusively as a passive vibration-suppression technology, a fundamental paradigm that the present study deliberately, rigorously, and systematically inverts.

Independently of surface modifications, the fundamental cross-sectional geometry of a bluff body exerts decisive and well-characterized influence on the scale of the aerodynamic forces it generates. The smooth circular cylinder (AR=1.0) serves as the canonical reference geometry throughout the VIV literature [16]. Elliptical cross-sections parametrized by the aspect ratio AR=b/a, where b is the streamwise semi-axis and a is the cross-flow semi-axis, represent a systematically controllable class of bluff bodies with markedly richer wake dynamics [17]. As the aspect ratio is reduced below unity while the frontal dimension a is held constant, the body presents an increasingly short and blunt streamwise profile to the oncoming flow. This enhanced bluffness compels the attached boundary layer to separate at a substantially earlier angular position from the forward stagnation point, elongating the free shear layer trajectory and enlarging the recirculation bubble in the near wake [18–20]. The expanded recirculation region is associated with a non-linear escalation of the mean pressure drag and, far more consequentially for VIV energy harvesting, a dramatic amplification of the fluctuating cross-flow lift force. For vertically-oriented bluff elliptical profiles, as adopted in the present study, this relationship is subject to a geometric limit: as the streamwise dimension becomes extremely short, the available boundary-layer development length diminishes and the shear-layer roll-up is constrained near the body surface, introducing a non-monotonic dependence of lift excitation on aspect ratio [21]. Furthermore, the Strouhal number does not vary monotonically with aspect ratio for the vertically-oriented geometry: the specific cross-sectional proportions govern the shear-layer formation frequency through a balance of separated flow length and vortex convection speed, which is examined in detail in the Section 3. These geometric effects make the present configuration particularly sensitive to the identification of an optimal aspect ratio, rather than simply minimising AR to maximise excitation [22,23].

Despite the substantial and parallel bodies of knowledge accumulated in biomimetic surface engineering and bluff-body aerodynamics, a critical and precisely defined gap persists at their intersection. Within the existing literature, grooved or textured surface modifications on cylindrical structures have predominantly been investigated as passive flow-control strategies for suppressing flow-induced vibrations and reducing structural loading [24–26]. Earlier studies on cactus-inspired and wavy cylinders primarily focused on the aerodynamic behaviour of these biomimetic surface modifications rather than on enhancing aerodynamic excitation for VIV energy harvesting applications [11,13,15]. Conversely, existing VIV energy harvesting research has concentrated almost exclusively on optimizing the mechanical and electromechanical transduction subsystem, including damping ratios, spring stiffnesses, piezoelectric coupling coefficients, and mass ratios, while treating the generator cylinder as a fixed, smooth, and circular geometric baseline [1,2,7,8,27]. No systematic computational optimization study has investigated whether biomimetic surface grooves, deliberately redesigned and parametrically optimized not for suppression but for amplification, can synergistically interact with a geometrically bluff elliptical base (AR<1.0) to maximize the available aerodynamic power for VIV energy harvesting. The potential synergistic interaction between optimized surface grooves and bluff elliptical cross-sections of varying aspect ratio, wherein the grooves perturb the separating shear layer and modulate vortex roll-up coherence while the elliptical bluffness governs the recirculation geometry and the available lift excitation magnitude, constitutes a relatively unexplored research area in passive aerodynamic excitation engineering. Critically, the interplay between groove spacing and boundary-layer thickness suggests that this synergy is not simply maximised by minimising the aspect ratio, but rather reaches an optimum at an intermediate bluffness level where the two geometric mechanisms are mutually reinforcing. The present study is explicitly motivated by, and systematically designed to address, this specific knowledge gap.

1.1 Objectives and Novelty of the Current Study

This study aims to explore how passive geometric modifications can amplify coherent wake forcing, and to quantify the resulting VIV energy harvesting potential for bluff bodies with tailored cross-sections. A fully coupled fluid-structure interaction model would inevitably introduce nonlinear feedback between cylinder motion and wake dynamics, complicating the isolation of purely geometric effects. To avoid this confounding influence, a fixed-cylinder computational framework is adopted, enabling the available aerodynamic excitation capacity to be evaluated strictly as a function of body geometry.

To achieve this objective, a three-stage computational framework is developed and executed within the laminar flow regime that corresponds directly to realistic ultra-low-speed hydrokinetic deployment conditions.

•   Stage I employs a Taguchi L9 orthogonal array at Re=200 to efficiently identify the optimal biomimetic groove wave number (N) and amplitude (Amp) that maximize the root-mean-square lift coefficient (Cℓ,rms) on a circular baseline (AR=1.0).

•   Stage II systematically couples the optimal groove profile with elliptical cross-sections of progressively decreasing aspect ratio (AR=1.0,0.75, and 0.50) at Re=200, quantifying how geometric bluffness interacts with the biomimetic surface modification to maximise aerodynamic excitation forces and near-wake fluctuation energy (WFE) concentration.

•   Stage III conducts a Reynolds-number sensitivity analysis at Re=100 and Re=150 for the two leading Grooved configurations (AR=0.75 and AR=0.50) and their smooth counterparts, establishing the robustness of the identified geometric optimum and the Reynolds-dependent activation of the boundary-layer saturation mechanism across the laminar deployment window.

The novelty of the present research comprises four distinct contributions:

•   Inverted biomimetic paradigm. Unlike previous studies that employ biomimetic surface grooves exclusively for vibration suppression, the present work deliberately redesigns and optimises groove geometry to amplify cross-flow aerodynamic excitation forces for energy harvesting purposes.

•   Geometric synergy concept. The coupled interaction between optimised surface grooves and bluff elliptical cross-sections of varying aspect ratio is systematically investigated, introducing the concept of geometric synergy to identify the optimal bluffness level at which the two geometric mechanisms are mutually reinforcing.

•   Multi-perspective energy harvesting quantification. POD modal analysis, WFE mapping, spectral (PSD) characterization, mean surface pressure distribution, and a simplified one-degree-of-freedom (1-DOF) structural projection model are integrated to provide a comprehensive, physics-informed assessment of the projected aerodynamic excitation potential for VIV energy harvesting across all evaluated configurations.

•   Reynolds-dependent geometric optimum. By extending the analysis to lower Reynolds numbers, the study demonstrates that the saturation of the groove contribution at AR=0.50 is a robust geometric feature of the laminar regime. In contrast, the absolute geometric optimum shifts from AR=0.50 at Re=100 to AR=0.75 at Re=200, indicating the onset of the boundary-layer saturation mechanism at Re=150.

2  Methodology

This section establishes the comprehensive computational and mathematical framework employed to evaluate the geometric amplification of aerodynamic excitation and the associated near-wake characteristics. To systematically isolate and quantify the impact of tailored geometric bluffness and biomimetic surface modifications without the confounding influence of structural feedback, a robust two-dimensional (2D) numerical methodology on fixed bodies is utilized throughout. The subsequent subsections detail: (i) the mathematical definition of the selected circular and elliptical geometries and the biomimetic groove parameterization; (ii) the governing fluid dynamics equations and the complete numerical setup within the laminar flow regime; (iii) the comprehensive independence studies and code validation procedures; and (iv) the theoretical foundations of both the Taguchi optimization method and Proper Orthogonal Decomposition (POD) modal analysis. This structured presentation is intended to ensure the full reliability, transparency, and reproducibility of all reported findings.

2.1 Geometric Modeling and Aerodynamic Formulations

The baseline circular and elliptical cross-sections are defined by the aspect ratio (AR):

AR=ba(1)

where a and b denote the full vertical (cross-flow) and horizontal (streamwise) axis lengths of the cross-section, respectively. Throughout this study, a is held strictly constant across all configurations, ensuring an identical frontal projection to the incoming flow. Accordingly, D=a serves as the reference length for all non-dimensional aerodynamic force and frequency calculations in the present study.

To achieve aspect ratios AR≤1.0, only b is systematically reduced. Three elliptical bases are evaluated, namely AR=1.0 (circular baseline), AR=0.75, and AR=0.50. Reducing b below a produces progressively shorter and more blunt streamwise bodies that impose increasingly severe adverse pressure gradients on the attached boundary layer, forcing earlier flow separation and broader near-wake recirculation regions. This deliberate geometric bluffness is the primary aerodynamic amplification mechanism exploited in Stage II.

The biomimetic grooved surface is generated by superimposing a cosinusoidal perturbation onto the base cross-section along the outward-facing surface normal. In a general vectorial form valid for all aspect ratios:

r(s)=r0(s)+a⋅AR⋅Amp⋅cos(2πNsL)n^(s)(2)

where r0(s) is the position vector on the smooth base ellipse, s is the arc-length coordinate measured from the leading stagnation point, L is the total perimeter of the base ellipse, N is the integer wave number, Amp is the non-dimensional groove amplitude expressed as a percentage of the streamwise dimension (b=a⋅AR), and n^(s) is the outward-facing unit normal vector. The arc-length parameterisation ensures physically uniform groove wavelength regardless of local curvature variations, and the proportional amplitude scaling with b prevents excessive structural distortion at lower aspect ratios.

For the circular baseline (AR=1.0, b=a), the normal coincides with the radial direction and Eq. (2) reduces to:

r(θ)=a2[1+2Amp⋅cos⁡(Nθ)](3)

where θ is the angular position from the leading stagnation point. Here the prefactor a/2=D/2 is the base radius of the circular cross-section, consistent with this full-dimension convention. Higher values of N produce more groove peaks per circumference, while larger Amp generates deeper surface undulations. The interaction between these two parameters governs the degree of shear layer perturbation, and their joint optimisation via the Taguchi method (Section 2.4.1) constitutes Stage I of the computational framework. Fig. 1 illustrates the baseline circular cross-section (AR=1.0) and the biomimetic grooved profile (N=24, Amp=5.0%). The sinusoidal surface undulation prescribed by Eq. (3) is clearly visible as a periodic radial perturbation superimposed on the circular baseline.

Figure 1: Schematic definition of the baseline circular (AR=1.0) and grooved (N=24, Amp=5.0%) geometries, illustrating the full axes a and b.

Additionally, the groove amplitude (Amp) is proportionally scaled with respect to the modified horizontal thickness b across different aspect ratios. This proportional scaling prevents excessive structural distortion of the grooved profiles at lower aspect ratios, where a fixed absolute amplitude would physically transform the surface undulations from the intended biomimetic cactus-inspired groove geometry into an extreme gear-tooth profile that would fundamentally alter the targeted boundary layer separation mechanisms.

The flow regime is characterized by the Reynolds number, defined using the constant frontal height D=a as the reference length:

Re=U∞Dν(4)

where U∞ is the uniform free-stream velocity and ν is the kinematic viscosity of the fluid. All simulations are conducted at this fixed Reynolds number, ensuring that the flow remains laminar and that any observed differences in aerodynamic behavior are solely attributable to geometric modifications rather than changes in flow regime.

The primary aerodynamic quantity of interest for VIV energy harvesting evaluation is the instantaneous cross-flow lift coefficient, Cℓ, defined as:

Cℓ=FL12ρU∞2D(5)

where FL is the instantaneous transient cross-flow lift force per unit span acting on the fixed bluff body, ρ is the fluid density, U∞ is the free-stream reference velocity, and D is the constant characteristic frontal height. The root-mean-square fluctuation of the lift coefficient, Cℓ,rms, quantifies the time-averaged magnitude of the alternating aerodynamic excitation force and serves as the primary optimization objective in this study:

Cℓ,rms=1T∫0T(Cℓ−Cℓ¯)2dt(6)

where Cℓ¯ denotes the time-mean lift coefficient and T is the integration interval encompassing a statistically sufficient number of fully developed, quasi-periodic vortex shedding cycles. A higher Cℓ,rms directly indicates a stronger alternating aerodynamic forcing, which in a coupled VIV system would drive larger oscillation amplitudes and correspondingly greater energy extraction capacity.

The dominant vortex shedding frequency fs is non-dimensionalized as the Strouhal number:

St=fsDU∞(7)

The Strouhal number provides a universal non-dimensional characterization of the vortex shedding periodicity across different geometric configurations and flow conditions. Changes in St with aspect ratio carry important structural implications: a lower St at a given Re translates to a lower dimensional shedding frequency fs, enabling structural lock-in to be achieved with softer mechanical systems or at lower flow velocities, both desirable properties for micro-scale VIV energy harvesters deployed in extremely slow currents. For the vertically-oriented elliptical geometry evaluated in this study, the Strouhal number does not vary monotonically with aspect ratio. The specific cross-sectional shape and the balance between shear-layer separation angle and vortex convection distance jointly determine the shedding frequency, producing a non-trivial dependence on AR that is characterised quantitatively in the Results section.

To evaluate the feasibility of utilizing the enhanced wake instabilities for practical energy extraction, a simplified one-degree-of-freedom (1-DOF) structural projection model is employed to estimate the available power coefficient (Cp,est):

Cp,est=Pest12ρU∞3D=Frms2/(2mζωn)12ρU∞3D(8)

where Pest is the projected extractable mechanical power, Frms is the RMS lift force derived directly from the fixed-body simulations (Frms=12ρU∞2D⋅Cℓ,rms), m is the effective oscillating system mass, ζ is the total structural-plus-electrical damping ratio, and ωn is the natural angular frequency of the theoretical coupled oscillator. In the present study, representative fixed values of the mass ratio m∗=m/md (where md is the displaced fluid mass) and the damping ratio ζ are adopted throughout, as detailed in Section 3.5.2. The 1-DOF projection model, originally rooted in the linear resonant response framework established by Bernitsas et al. [1] and widely adopted in the VIV energy harvesting literature [28–31], provides a standardized and physically meaningful metric that enables direct comparison of the geometric amplification potential across different configurations, independent of the specific structural or electromechanical transduction system ultimately employed. The present 1-DOF model provides a relative measure of aerodynamic excitation across configurations rather than an absolute prediction of harvested power. Since the structural parameters are held fixed, Pest is interpreted throughout as a comparative fixed-body aerodynamic forcing metric. It quantifies the relative excitation available to a hypothetical coupled oscillator and reflects comparative trends across configurations. It is not a prediction of the actual harvested power, oscillation amplitude, or lock-in behaviour, which depend on the specific device design and would require fully coupled fluid-structure interaction simulations.

2.2 Numerical Setup and Flow Configuration

All transient numerical simulations are performed using a two-dimensional, pressure-based, laminar flow formulation within the ANSYS Fluent commercial CFD solver environment. The simulation framework is deliberately two-dimensional, which accurately captures the dominant planar Kármán vortex street dynamics and the associated periodic aerodynamic force fluctuations at Re = 200 while affording the computational efficiency necessary for the parametric optimization and comparative studies constituting the core of this investigation. Although secondary three-dimensional wake instabilities (Mode A transition) begin to emerge in circular cylinder wakes at Re around 190 [32], numerical studies conducted specifically at Re = 200 have demonstrated that two-dimensional simulations reproduce the aerodynamic force coefficients and Strouhal number with negligible deviation from their three-dimensional counterparts [33–35], confirming that the primary energy-carrying Kármán vortex street and its associated cross-flow forcing remain fundamentally two-dimensional at this Reynolds number. Although these three-dimensional effects, principally spanwise vortex dislocations, would drain a fraction of the spanwise-coherent fluctuation energy and somewhat reduce the absolute force magnitudes, they act on the same primary Kármán mode shared by all six geometries; the comparative ranking and amplification factors that constitute the central conclusions of this study are expected to remain qualitatively similar, although some quantitative changes in amplification factors and relative ranking cannot be excluded without dedicated three-dimensional simulations. The Reynolds-sensitivity analysis at Re = 100 and Re = 150 lies, moreover, entirely within the strictly two-dimensional laminar regime. Since the present study targets comparative aerodynamic excitation across geometries rather than absolute force prediction, the two-dimensional framework is well suited to the intended geometric optimisation, with the reported absolute force levels interpreted as upper-bound estimates.

The two-dimensional, incompressible, unsteady laminar flow is governed by the continuity equation and the Navier-Stokes momentum equations:

∇⋅u=0(9)

∂u∂t+(u⋅∇)u=−1ρ∇p+ν∇2u(10)

where u is the velocity vector field comprising the streamwise (u) and cross-flow (v) components, t denotes physical time, p is the pressure field, ρ is the fluid density, and ν is the kinematic viscosity.

To ensure robust and accelerated convergence for the highly unsteady wake dynamics, the Pressure-Implicit with Splitting of Operators (PISO) pressure-velocity coupling algorithm is employed, utilizing both skewness and neighbor corrections to maintain solution stability on the complex unstructured mesh topology near the grooved cylinder surfaces. Spatial discretization of velocity gradients is performed using the Green-Gauss Node Based method. A second-order scheme is applied for pressure interpolation, while the momentum equations are discretized using a second-order upwind scheme to minimize artificial numerical diffusion across the computational grid [36]. Temporal advancement is governed by a second-order implicit (Crank-Nicolson blended) formulation, which precisely captures the time-dependent, quasi-periodic vortex shedding dynamics and the resulting aerodynamic force fluctuations with minimal phase error [37].

The computational domain is a two-dimensional rectangular channel surrounding the fixed cylinder model. The domain boundaries are parameterized using three key distances: the upstream length A from the velocity inlet to the cylinder center, the downstream length B from the cylinder center to the pressure outlet, and the total lateral height C of the domain, as illustrated in Fig. 2. A uniform, steady velocity inlet condition (U∞=const., v=0) is applied at the upstream boundary. A zero gauge pressure outlet condition is imposed at the downstream boundary, permitting the fully developed vortex wake to exit the domain without artificial reflection. Symmetry boundary conditions are applied at the upper and lower lateral boundaries to simulate an unconfined far-field environment. A no-slip, no-penetration wall boundary condition (u=v=0) is applied across the entire cylinder surface. The cylinder is stationary and fixed throughout all simulations.

Figure 2: Schematic of the computational domain illustrating the upstream distance (A), downstream wake region (B), lateral height (C), and applied boundary conditions.

The computational mesh consists of a structured quadrilateral cell topology throughout the far-field domain, transitioning to a carefully designed multi-layer boundary layer mesh in the near-wall region. The boundary layer mesh is constructed with geometric stretching, with the first computational cell layer set at a thickness of 6×10−4 m from the wall, directly resolving the viscous boundary layer characteristic of low-Reynolds-number laminar flows without requiring any wall function treatment.

Once the transient simulations reach a statistically steady, fully developed periodic state, confirmed by the stabilization of the lift coefficient time history to a consistent quasi-sinusoidal signal, the initial start-up phase (the first 50% of the total simulation time) is strictly discarded to prevent any transient contamination. For the fully developed regime, instantaneous velocity field snapshots are exported from the solver at intervals of 50 computational time steps (Δtsample=0.5 s). These snapshots comprise the streamwise and cross-flow velocity components at all grid nodes. This sampling rate is deliberately chosen to yield approximately 10 snapshots per shedding cycle, providing more than twice the minimum sampling frequency required by the Nyquist–Shannon sampling theorem [38] for the dominant Strouhal shedding frequency while maintaining efficient computational storage and memory allocation. The extracted temporal window encompasses 20 consecutive fully developed vortex shedding cycles, yielding a matrix of 200 snapshots. This dataset constitutes the velocity field input required for the subsequent Proper Orthogonal Decomposition (POD) modal analysis and wake fluctuation energy (WFE) spatial mapping.

2.3 Independence Studies and Model Validation

To guarantee the numerical accuracy and physical fidelity of all computed results, and to eliminate any contamination from spatial or temporal discretization errors, a comprehensive three-part independence study is conducted on the smooth circular cylinder at Re=200 prior to performing any geometric modification simulations. The three independence parameters evaluated are: (i) the computational domain extent, (ii) the mesh spatial resolution, and (iii) the physical time step size. Each parameter is varied systematically while keeping the other two fixed at their refined levels.

2.3.1 Domain Size Independence

The spatial independence study is performed first to eliminate unphysical boundary interference effects that arise when the domain boundaries are placed too close to the cylinder. Four domain configurations of progressively increasing extent are evaluated, as detailed in Table 1. The upstream distance A is chosen to ensure that the inlet boundary condition is imposed in a region of genuinely undisturbed, uniform flow, while the downstream distance B must be sufficient to allow the periodic vortex street to fully develop and dissipate before reaching the pressure outlet.

The aerodynamic coefficients, namely the mean drag coefficient (Cd¯) and RMS lift (Cℓ,rms), are monitored and compared across the four domain sizes, as shown in Fig. 3. The results demonstrate complete stabilization of both quantities at Domain 3 (30D×90D×60D), beyond which further domain expansion to Domain 4 produces changes of less than 0.1% in both Cd¯ and Cℓ,rms. Domain 3 yields a solid blockage ratio of approximately 1.67% (D/C=1/60), which is enough for laminar bluff body simulations and effectively eliminates artificial far-field confinement effects [39]. Domain 3 is therefore selected as the standard computational domain for all subsequent simulations.

Figure 3: Domain size independence study.

2.3.2 Grid Resolution Independence

Following the domain selection, a rigorous grid resolution independence study is conducted to determine the minimum spatial discretization required for accurate resolution of the near-wall boundary layer and the shear layer dynamics. Four structured mesh configurations of progressively increasing density are generated, differing in the number of computational nodes distributed around the cylinder circumference, as summarized in Table 2. In all mesh configurations, the first layer thickness and the boundary layer expansion ratio (1.1) are maintained constant, ensuring that only the circumferential node count governs the differences between refinement levels.

The grid convergence behaviour of Cd¯ and Cℓ,rms across the four mesh levels is presented in Fig. 4. The transient aerodynamic force convergence analysis demonstrates that both Cd¯ and Cℓ,rms become fully grid-independent beyond Mesh 3 (240 nodes around circumference), with differences between Mesh 3 and the finest Mesh 4 remaining below 0.3%. Mesh 1 substantially underestimates Cℓ,rms due to insufficient resolution of the shear layer curvature near the separation point. The substantially greater computational cost of Mesh 4 relative to Mesh 3, combined with its negligible improvement in predicted aerodynamic coefficients, justifies the selection of Mesh 3 as the standard computational mesh for all production simulations.

Figure 4: Grid size independence study.

The grid-independence study established a converged baseline resolution of 240 circumferential nodes (Table 2). For the grooved and elliptical geometries, the mesh is scaled to preserve the baseline cell size while adequately resolving the grooves, with identical wall-normal first-cell thickness across all cases. Consequently, all configurations inherit the grid-independent baseline resolution, with finer effective resolution in the groove regions. The final computational mesh employed in the simulations is shown in Fig. 5, comprising a structured quadrilateral far-field grid and a refined near-wall boundary layer mesh around the baseline and grooved cylinders.

Figure 5: Computational mesh topology illustrating the structured far-field quadrilateral grid and the highly refined near-wall boundary layer mesh around the baseline and grooved circular cylinder.

2.3.3 Temporal Discretization Independence

The temporal resolution study evaluates the sensitivity of the computed transient aerodynamic forces to the physical time step size Δt. A sufficiently small time step (TS) is required to accurately resolve the instantaneous dynamics of vortex formation and shedding, satisfy the Courant-Friedrichs-Lewy (CFL) stability criterion within the highly unsteady near-wake region, and ensure that the spectral content of the lift force signal is captured without aliasing or numerical damping. Three time step sizes are evaluated, as listed in Table 3.

The independence test results, as illustrated in Fig. 6, confirm that the transient force fluctuations, including their amplitude, phase, and spectral content, are accurately and consistently captured at a time step of Δt=0.01 s (TS 2). This temporal resolution ensures that the maximum CFL number remains below 1.0 throughout the critical near-wake and boundary layer regions. The coarser time step TS 1 (Δt=0.02 s) introduces measurable numerical damping that artificially reduces the peak-to-peak lift amplitude by approximately 4%, while TS 3 (Δt=0.005 s) produces negligible changes relative to TS 2 while incurring approximately twice the computational cost.

Figure 6: Time step size independence study.

Accordingly, the combination of Domain 3, Mesh 3, and Time Step Size 2 is selected as the standard computational configuration for all transient simulations. Each simulation is executed for a minimum of 20,000 computational time steps (200 s of physical time), 34–44 complete vortex shedding cycles depending on the specific Strouhal frequency of the simulated configuration.

2.3.4 Code Validation

The physical and numerical parameters employed throughout the study are summarised in Table 4 for completeness and reproducibility. The validation of the adopted numerical framework is then performed using the baseline smooth circular cylinder case.

To validate the numerical framework and solver settings, the baseline smooth circular cylinder is simulated at Re = 200 and the resulting aerodynamic force coefficients and Strouhal number are compared against well-established reference data, comprising the direct numerical simulation of Henderson [40], the experimental measurements of Williamson [41], the lattice Boltzmann simulation of Hu et al. [34], and the finite-volume study of Rajani et al. [35].

The present simulation predicts a mean drag coefficient of Cd¯=1.333 in good agreement with Henderson [40] and Rajani et al. [35], an RMS lift coefficient Cℓ,rms=0.4765 within the range reported by Hu et al. (0.499) and Rajani et al. (0.4244), and a Strouhal number St=0.1953 that matches the reference values to within their reported uncertainty bounds, as shown in Table 5. Notably, the two-dimensional Strouhal number agrees closely with the experimentally measured value of Williamson [41], whose physical flow is inherently three-dimensional, indicating that the primary Kármán shedding characteristics are captured reasonably well by the adopted two-dimensional framework at Re = 200. The framework is therefore reliable for the subsequent grooved and elliptical configurations, and all aerodynamic metrics are presented and interpreted as comparative trends across geometries within this two-dimensional laminar setting.

2.4 Optimization Framework and Modal Analysis

2.4.1 Taguchi L9 Orthogonal Array Optimization

To systematically and efficiently identify the biomimetic groove configuration that yields the maximum aerodynamic excitation without incurring the prohibitive computational cost of a full factorial parametric study, the Taguchi robust design optimization method is employed in the first stage of the computational framework. The Taguchi method utilizes orthogonal arrays to strategically sample the design space in a balanced and resolution-maximizing manner [42,43]. An L9 orthogonal array is implemented in this study, evaluating three levels of each of the two groove parameters (N and Amp) in nine strategically selected simulation runs. The L9 array provides full resolution of both main effects while requiring the minimum number of simulation runs needed to achieve this resolution. The specific parameter levels evaluated are: wave number N∈{20,24,28} and amplitude Amp∈{3%,4%,5%} of the reference length D. The lower bound of the N range is set at 20 to ensure a minimum groove wavelength sufficient for meaningful shear layer interaction, while the upper bound of 28 represents the practical limit beyond which the inter-groove angular spacing approaches the length scale of the near-wall viscous sublayer at the present laminar Reynolds number. The amplitude range balances the requirement for a geometrically significant surface perturbation (Amp≥3%) against the risk of excessive structural distortion (Amp>5%), which would compromise the biomimetic character of the grooved profile.

Because the primary objective for VIV energy harvesting is to maximize the alternating cross-flow aerodynamic force, the “Larger-the-Better” Signal-to-Noise (S/N) ratio is selected as the Taguchi objective:

S/N=−10log10[1n∑i=1n1yi2](11)

where yi corresponds to the root-mean-square lift coefficient (Cℓ,rms) obtained from each individual simulation run i in the L9 array, and n is the number of observations per run. Since the laminar flow at fixed Reynolds number is deterministic and reproducible, each configuration was simulated once, and n corresponds to the single converged Cℓ,rms obtained from a statistically stationary lift signal. The L9 array therefore serves as an efficient structured design of experiments for the main effects of N and Amp, with the S/N ratio acting as a monotonic objective transform of Cℓ,rms rather than a noise-robustness measure. The associated numerical uncertainty is bounded by the independence studies of Section 2.3, below approximately 0.3% in Cℓ,rms, with correspondingly small uncertainties in St and Pest. The main effects of each parameter level are subsequently assessed by averaging the S/N ratios across all runs sharing a given parameter level, enabling the optimal combination to be identified as the parameter levels that individually maximize the average S/N response.

2.4.2 Proper Orthogonal Decomposition (POD) Modal Analysis

Following the identification and coupling of the optimal groove geometry with the elliptical base configurations, Proper Orthogonal Decomposition (POD) is applied to the time-resolved near-wake velocity field data to elucidate the underlying physical mechanisms driving the aerodynamic force amplification and to systematically quantify the spatial coherence and energetic organization of the vortex shedding structures [44,45]. POD is a mathematically optimal technique for extracting the most energetically dominant spatial flow structures from an ensemble of velocity field snapshots: it decomposes the highly unsteady near-wake into a ranked series of orthogonal spatial modes ordered strictly by their wake fluctuation energy content, ensuring that the first mode always captures the greatest possible fraction of the total fluctuation energy.

The snapshot-based POD method is implemented [46]. The fluctuating velocity field u′(x,t) is expressed as a linear superposition of M time-independent orthogonal spatial modes φk(x) and their corresponding time-varying scalar coefficients ak(t):

u′(x,t)=∑k=1Mak(t)φk(x)(12)

The spatial modes φk are computed as the eigenvectors of the two-point spatial cross-correlation matrix constructed from the snapshot ensemble, while the time coefficients ak(t) are obtained by projecting each snapshot onto the corresponding spatial mode. The eigenvalue λk associated with each mode directly represents the energy of wake fluctuation contributed by that mode, and the normalized fractional energy content of each mode is given by λk/∑λ. The cumulative energy fraction provides a rigorous quantitative measure of the degree of wake coherence: for a highly organized, synchronized vortex shedding process, as expected and desired for sustained VIV energy harvesting, the first two spatial modes together capture the overwhelming majority (typically >95%) of the total fluctuation energy, confirming the existence of a dominant, low-dimensional, periodic wake structure.

The physical significance of POD in the context of VIV energy harvesting is profound: a wake dominated by its first two spatial modes is a wake in which virtually all of the available fluctuation energy is organized into a single coherent, periodic oscillation at the Strouhal frequency. This high modal coherence implies that the aerodynamic lift force on the cylinder is similarly dominated by a single, stable periodic component, representing the forcing quality required to drive sustained lock-in oscillations in a coupled structural resonator. The POD modal energy distribution therefore serves as a critical diagnostic for assessing not merely the magnitude of the aerodynamic excitation, but also its regularity and temporal stability, both essential qualities for reliable, long-duration energy harvesting in practical deployment.

3  Results and Discussion

The results are organised into six thematic subsections following the logical progression of the computational framework. The Taguchi L9 optimisation outcome and optimal groove selection are presented first in Section 3.1, followed by the aerodynamic force and spectral response characterisation in Section 3.2. The modal coherence and structural stability diagnostics are presented in Section 3.3, and near-wake topology and instantaneous dynamics are examined in Section 3.4. The energy-oriented wake assessment is presented in Section 3.5, where the projected energy harvesting amplification is quantified across all six configurations and the optimum geometry is identified. The study culminates in the Reynolds-number sensitivity analysis in Section 3.6, which establishes the robustness and Reynolds dependence of the identified geometric optimum across the laminar regime.

3.1 Parametric Optimisation of Biomimetic Grooves

The Taguchi L9 simulation matrix and the resulting aerodynamic metrics are presented in Table 6. The signal-to-noise ratio analysis identifies groove amplitude (Amp) as the dominant control parameter. Averaging the S/N ratios across all runs at each amplitude level confirms a monotonic and physically consistent improvement with increasing groove depth, whereas the wave number N exerts a negligible main effect, with the level-averaged S/N ratios for N=20, 24, and 28 differing by less than 0.02 dB, indicating that the influence of N is substantially weaker than that of Amp over the parameter range considered. This outcome agrees with previous laminar studies showing that shear-layer instability is governed primarily by perturbation amplitude rather than spatial frequency [13].

The physical origin of this asymmetry lies in the laminar near-wall dynamics at Re=200. A deeper groove imposes a stronger localised adverse pressure gradient at each ridge, accelerating shear-layer detachment and intensifying the subsequent vortex roll-up regardless of circumferential density. The spatial frequency N, by contrast, governs the inter-groove angular spacing; at N=24 this spacing is approximately 15∘, which is well matched to the coherent near-wall flow structures at this Reynolds number. Further increasing N to 28 tightens the spacing to 12.9∘, at which point the individual groove ridges begin to act collectively as a near-wall roughness layer rather than as discrete shear-layer tripping elements, diminishing their aerodynamic effectiveness [13,47].

With the dominant parameter fixed at Amp=5%, three candidate configurations remain, namely Run 03 (N=20), Run 06 (N=24), and Run 09 (N=28). Runs 06 and 09 yield virtually identical Cℓ,rms values (0.5768 vs. 0.5766), with a difference of less than 0.04%, which is negligible for practical purposes. The selection is therefore governed by physical consistency rather than numerical ranking. At N=24 the inter-groove spacing preserves the intended biomimetic tripping character, whereas N=28 risks shifting below the effective groove-to-boundary-layer interaction range established for cactus-inspired geometries [13,47]. Run 06 (N=24, Amp=5%) is accordingly selected as the optimal groove configuration for Stage II.

Run 06 achieves Cℓ,rms=0.5768 and Pest=0.3873 W/m, representing a 21% increase in Cℓ,rms and a 12.5% downward shift in Strouhal number (from St=0.1953 to St=0.1709) relative to the smooth baseline. This Strouhal reduction is structurally advantageous for low-speed lock-in, as it lowers the natural frequency required for synchronisation and broadens the operating range of a micro-harvester in ultra-low-velocity environments [1,8].

The groove parameters (N and Amp) were optimised on the circular baseline at Re = 200 and subsequently transferred to the elliptical and lower-Reynolds-number configurations without re-optimisation. Given the weak sensitivity to N identified by the Taguchi analysis and the fact that Amp=5% was the largest amplitude considered and yielded the highest aerodynamic performance within the examined range, the selected groove geometry was retained for all subsequent cases.

3.2 Aerodynamic Force Response and Spectral Characteristics

3.2.1 Force Coefficient Comparison

Fig. 7 presents the mean drag coefficient (Cd¯) and RMS lift coefficient (Cℓ,rms) for all six Stage II configurations, grouped by aspect ratio with smooth baseline and optimally grooved cases displayed side by side.

Figure 7: Mean drag (Cd¯) and RMS lift (Cℓ,rms) coefficients for the smooth and optimally grooved configurations across varying aspect at Re = 200.

Examining the smooth baseline series first, increasing geometric bluffness from AR=1.0 to AR=0.50 raises Cd¯ monotonically from 1.333 to 1.891, confirming the expected non-linear escalation of pressure drag as the streamwise body dimension is shortened and an increasingly severe adverse pressure gradient is imposed on the attached boundary layer. The RMS lift coefficient follows a different trend: Cℓ,rms increases from 0.4765 at AR=1.0 to 0.6314 at AR=0.75, but then falls slightly to 0.6126 at AR=0.50. This non-monotonic behaviour is a characteristic of the vertically-oriented elliptical geometry adopted in this study. At AR=0.50, the streamwise body is so compressed that the boundary layer barely develops before separating, limiting the effective shear-layer roll-up length and partially constraining the alternating lift magnitude despite the increased drag bluffness. This physical mechanism is examined in detail through the wake topology analysis in Section 3.4.

The introduction of biomimetic grooves (N=24, Amp=5%) consistently elevates both Cd¯ and Cℓ,rms relative to the corresponding smooth baseline at every aspect ratio, confirming that the optimised surface modification successfully inverts the conventional vibration-suppression role of such grooves. Three aspects of this groove contribution are particularly significant. First, at AR=1.0 the groove raises Cℓ,rms by 21.1% (from 0.4765 to 0.5768), representing the largest relative increment across all three aspect ratios. Second, at AR=0.75 the groove produces a 13.0% increment (from 0.6314 to 0.7134), yielding the highest absolute Cℓ,rms of any evaluated configuration. Third, at AR=0.50 the groove contribution diminishes sharply to only 7.0% (from 0.6126 to 0.6555). This monotonically decreasing groove effectiveness with increasing bluffness reveals a physical saturation: as the elliptical body becomes progressively shorter in the streamwise direction, the boundary layer available for groove interaction thins and shortens, reducing the capacity of the ridge-induced micro-vortices to coherently perturb the separating shear layer [13,47]. This behaviour suggests that the effectiveness of groove-induced perturbations is constrained by the available boundary-layer development length, which decreases with increasing geometric bluffness. The critical implication is that the highest absolute Cℓ,rms=0.7134 is achieved not at the maximum bluffness level but at the intermediate Grooved AR=0.75 configuration, identifying this geometry as the optimal configuration within the investigated parameter space.

To provide surface pressure evidence supporting the observed force coefficient trends, Fig. 8 presents the time-averaged mean pressure coefficient (Cp¯) distribution around the cylinder circumference for four representative configurations. The inset schematic defines the angular coordinate convention, with θ=0∘ at the leading stagnation point and θ=180∘ at the base region. Faint background traces represent instantaneous groove-induced pressure perturbations, while bold markers indicate the circumferentially averaged mean values.

Figure 8: Time-averaged mean pressure coefficient (Cp¯) distribution for baseline AR=1.0, grooved AR=1.0, grooved AR=0.75, and grooved AR=0.50 at Re = 200. Faint traces indicate instantaneous groove-induced perturbations. Inset defines the angular coordinate convention.

Baseline AR=1.0 follows the canonical laminar profile with Cp¯=1.0 at the stagnation point and base pressure of approximately −0.7 at θ=180∘. The Grooved AR=1.0 configuration exhibits a moderately deeper suction envelope with organized periodic pressure perturbations associated with coherent shear-layer tripping by the groove ridges on the separating shear layer. At Grooved AR=0.75, the suction amplifies substantially, with Cp¯ in the range of −1.6 to −2.0 across the leeward surface, and the groove-induced oscillations remain spatially organised and periodic, confirming coherent groove-shear-layer interaction at this bluffness level. This large and coherent pressure asymmetry between the upper and lower surfaces directly generates the highest Cℓ,rms=0.7134 of any evaluated configuration.

In contrast, Grooved AR=0.50 produces the most extreme suction peak, consistent with its highest Cd¯=1.987, yet the groove-induced perturbations become visibly irregular and incoherent in the separation regions. At this extreme bluffness, the boundary layer traverses the groove ridges too rapidly for the biomimetic tripping mechanism to establish organised micro-vortex structures, producing incoherent pressure fluctuations that less efficiently convert the high suction levels into coherent alternating lift. The Cp¯ distribution therefore provides an independent line of evidence confirming the Grooved AR=0.75 configuration as the aerodynamically optimal configuration among the evaluated cases.

3.2.2 Spectral Analysis and Strouhal Characteristics

Fig. 9 presents the Welch power spectral density (PSD) of the lift coefficient signal on a logarithmic scale for all six configurations, arranged as three panels each directly comparing the baseline (solid line) and grooved (dashed line) spectra at a common aspect ratio. All spectra are dominated by a single sharp peak at the Kármán shedding frequency, confirming the highly periodic, narrow-band character of the laminar wake across all evaluated geometries.

Figure 9: Power spectral density (PSD) of lift fluctuations comparing baseline and grooved configurations at AR=1.0, 0.75, and 0.50 at Re = 200.

At AR=1.0 (upper panel), the groove induces a 12.5% downward shift in the dominant spectral peak from St=0.1953 to St=0.1709. This is the largest Strouhal reduction observed across all configurations. The groove ridges promote early shear-layer detachment, which effectively shortens the vortex formation length and reduces the convection speed of the shed vortices relative to the free stream, lowering the dimensional shedding frequency. From a practical harvesting perspective, this frequency reduction translates directly into a lower natural frequency requirement for structural lock-in, which is advantageous for the design of soft, low-stiffness resonators suited to ultra-low-speed hydrokinetic deployment [7,8].

At AR=0.75 (middle panel), the groove produces a Strouhal normalisation from St=0.2197 (smooth) to St=0.1953. This restoration of the shedding frequency to the circular cylinder reference value carries a significant practical advantage: a VIV micro-harvester designed and tuned for the standard circular cylinder geometry could in principle be deployed with the Grooved AR=0.75 body without any structural retuning of the spring-mass resonator, making this configuration especially attractive for scalable device implementation [1].

At AR=0.50 (lower panel), the Strouhal number remains at St=0.1953 in both the baseline and grooved cases, indicating that the shedding frequency at this extreme bluffness level is dominated by the geometry of the separated shear layers and is insensitive to the groove perturbation. In all panels, a secondary spectral peak is visible at approximately St=0.59 in the AR=0.50 and AR=0.75 configurations, corresponding to the second harmonic of the fundamental Kármán frequency. This feature is a well-established characteristic of laminar bluff-body wakes at Re=200 and is consistent with the self-interaction of the periodic vortex street in the near-wake region [40].

3.3 Modal Coherence and Structural Stability

3.3.1 POD Modal Energy Distribution

Fig. 10 presents the POD modal energy fraction for the first four modes and the cumulative energy fraction for the same four representative configurations. The data provide a quantitative measure of wake coherence that informs the forcing quality of each configuration relevant to sustained VIV energy harvesting. For Baseline AR=1.0, Mode 1 and Mode 2 together capture 96.5% of the total fluctuation energy (54.6+41.9%), with Modes 3 and 4 accounting for only 1.4% and 0.9% respectively. This near-perfect two-mode dominance confirms the textbook low-dimensional character of the laminar circular cylinder wake at Re=200, in which virtually all aerodynamic fluctuation energy is organised into the conjugate Kármán mode pair. The Grooved AR=1.0 configuration maintains essentially identical modal coherence, with the first two modes capturing 96.2% (54.7+41.5%). This near-invariance confirms that the groove modification at AR=1.0 intensifies the aerodynamic excitation without disrupting the fundamental low-dimensional periodicity of the wake, consistent with the modest but consistent Cℓ,rms increment observed in Fig. 7.

Figure 10: Proper orthogonal decomposition (POD) modal and cumulative energy distributions for the first four modes Re = 200.

At the optimum Grooved AR=0.75, the first two modes together capture 96.2% (54.5+41.7%) of the total energy, demonstrating that the geometric synergy between the groove and the intermediate bluffness level preserves near-perfect modal coherence while simultaneously achieving the highest Cℓ,rms of any evaluated configuration. From an energy harvesting perspective, this result is highly significant: the amplified aerodynamic forcing of the Grooved AR=0.75 configuration is entirely organised into a single dominant periodic oscillation at the Strouhal frequency. A mechanical resonator coupled to this wake would therefore be expected to experience a reliable, narrow-band harmonic forcing favourable for energy extraction through the lock-in mechanism [5].

A marked departure from this high-coherence regime occurs at Grooved AR=0.50. The combined energy fraction of Modes 1 and 2 drops to 88.1% (49.1+39.0%), accompanied by a sharp anomalous increase in the energy scattered into higher-order structures: Modes 3 and 4 each carry 3.9% of the total energy, more than double their contribution in all other configurations. This spectral leakage into higher modes demonstrates that the extreme adverse pressure gradients associated with the AR=0.50 bluffness level overwhelm the organising influence of the biomimetic grooves, forcing the wake into a higher-dimensional, more complex dynamical state. As a consequence, the aerodynamic excitation at AR=0.50 loses the optimal narrow-band quality that is essential for efficient VIV energy conversion, fully elucidating why Cℓ,rms does not increase monotonically with further bluffness and confirming that the AR=0.75 configuration represents the optimum case in the present study.

3.3.2 POD Spatial Mode Structures

Fig. 11 presents the spatial contours of POD Modes 1 through 4 for the four evaluated configurations, providing a spatially resolved representation of the physical structures that carry the modal energy identified in the preceding subsection. The dominant vortex-street wavelength λ, annotated on the Mode 1 panels, is measured as the streamwise distance between two consecutive same-sign vorticity peaks. This wavelength is related to the shedding frequency through λ=Uc/fs, where Uc is the mean convection velocity of the vortical structures. Introducing the Strouhal number definition yields λ/a=(Uc/U∞)/St. For laminar cylinder wakes at Re=200, the convection velocity remains of the same order as the the free-stream velocity (Uc in the range of 0.7–0.9U∞) [48], so that λ scales approximately inversely with St. Accordingly, the wavelength variations across configurations are interpreted primarily in terms of relative changes in St.

Figure 11: Spatial contours of the first four POD modes for baseline AR=1.0, grooved AR=1.0, grooved AR=0.75, and grooved AR=0.50 Re = 200. The dominant vortex-street wavelength λ is annotated in the mode 1 panels.

For Baseline AR=1.0, Modes 1 and 2 display the canonical alternating positive-negative vorticity patterns of the Kármán street, propagating coherently downstream with a dominant wavelength (λ) of approximately 4.10D that is directly related to the Strouhal shedding frequency through the convection velocity. Modes 3 and 4 exhibit low-amplitude spatial structures confined to the immediate near-wake, confirming their role as higher-harmonic corrections with negligible energetic significance.

The Grooved AR=1.0 configuration produces a noticeably longer dominant spatial wavelength of λ around 4.68D in Modes 1 and 2, fully consistent with the lower Strouhal number (St=0.1709) induced by the groove modification. A lower shedding frequency at a given convection velocity corresponds directly to a larger spatial period between consecutive vortex cores, and the POD spatial structure provides independent geometric confirmation of the spectral shift observed in Fig. 9.

At Grooved AR=0.75, the dominant wavelength returns to approximately 4.10D, consistent with the Strouhal normalisation to St=0.1953 documented in Section 3.2. The Mode 1 and 2 structures at AR=0.75 display intensified near-field vorticity relative to both AR=1.0 cases, reflecting the amplified vortex roll-up associated with the shorter recirculation zone. Critically, the anti-symmetric Kármán pattern propagates coherently throughout the downstream domain with no visible fragmentation, confirming the sustained structural integrity of the wake at this optimal bluffness level.

In sharp contrast, the Grooved AR=0.50 configuration reveals a progressive spatial breakdown in Modes 1 and 2: while the anti-symmetric pattern is established in the near-wake, the structures rapidly break down into disjointed, irregular clusters as they convect downstream, losing spatial coherence well within the computational domain. More strikingly, Modes 3 and 4 at AR=0.50 exhibit massively extended spatial footprints reaching far into the downstream wake, directly visualising the incoherent multi-scale vortical structures that absorb the anomalous higher-mode energy fraction identified in Fig. 10. This spatial fragmentation confirms that the aerodynamic energy at AR=0.50 is partially dispersed into structures that cannot coherently excite a 1-DOF structural resonator, providing a physics-based explanation for the energy harvesting performance ceiling at this bluffness level.

3.4 Mean Wake Topology and Instantaneous Dynamics

3.4.1 Mean Recirculation Zone

Fig. 12 presents the time-averaged streamwise velocity contours for four configurations that span the full parameter space of Stage II: Baseline AR=1.0, Grooved AR=1.0, Grooved AR=0.75, and Grooved AR=0.50. The recirculation zone length Lf is annotated on each panel as the downstream distance from the cylinder base to the point where the time-mean streamwise velocity first recovers to zero along the wake centreline.

Figure 12: Mean wake topology and recirculation zone lengths (Lf) for the evaluated geometric configurations Re = 200.

For the smooth circular cylinder, Lf=0.87D, consistent with the established literature for laminar circular cylinder wakes at Re=200 [39]. The introduction of biomimetic grooves at AR=1.0 marginally contracts the recirculation zone to Lf=0.85D. This small but consistent reduction reflects the earlier shear-layer detachment promoted by the groove ridges, which slightly shortens the vortex formation region without fundamentally altering the global wake topology. This observation is consistent with prior findings for wavy cylinders in the laminar regime, in which surface undulations promote earlier separation through localised adverse pressure gradients without dramatically changing the time-averaged wake structure [13].

A fundamentally different topological regime emerges when the elliptical bluffness is introduced. At Grooved AR=0.75, the combined effect of the shortened streamwise body and the groove perturbations drives a drastic contraction of the recirculation zone to Lf=0.61D, representing a 30% reduction relative to the smooth circular baseline. This compact recirculation bubble implies that coherent vortex roll-up is completed significantly closer to the cylinder base. The vortices are therefore deposited into a high-intensity, near-body region that generates a correspondingly large alternating low-pressure suction force, directly driving the elevated Cℓ,rms=0.7134 of this configuration. At Grooved AR=0.50, the recirculation zone further contracts to Lf=0.56D, accompanied by a pronounced lateral broadening of the zero-velocity contour into a wide, flattened stagnation region. This topological change signals the onset of the saturation regime: the shear layers now detach so abruptly and symmetrically that the organised vortex roll-up process is partially disrupted, which accounts for the reduced Cℓ,rms relative to Grooved AR=0.75 despite the shorter Lf.

3.4.2 Instantaneous Vorticity Fields

Fig. 13 presents the instantaneous normalised vorticity contours for the same four configurations at a representative phase of the statistically steady shedding cycle. The qualitative evolution of the instantaneous wake structure across the four configurations provides direct visual evidence of the physical mechanisms identified in the force and topology data.

Figure 13: Instantaneous spanwise vorticity field at Re = 200, normalised as ωz∗=ωzD/U∞.

For Baseline AR=1.0, the wake consists of a well-organised Kármán vortex street in which compact, alternating positive and negative vorticity cores convect smoothly downstream at regular spatial intervals, maintaining their structural integrity well into the far wake. The Grooved AR=1.0 configuration preserves this globally organised character but produces slightly more intensified and spatially compact vortex cores in the near-wake region, consistent with the groove-induced acceleration of the roll-up process and the modest Cℓ,rms increment of 21.1%.

At the optimum Grooved AR=0.75, the instantaneous vorticity field transforms markedly. The shear layers separate from the wider, shorter elliptical body and roll up into distinctly swollen, high-intensity vortex cores that form at a reduced downstream distance from the cylinder base. Despite the intensified vorticity magnitude, the shedding pattern retains a highly organised, anti-symmetric alternating character, confirming that the AR=0.75 bluffness level lies within the regime of coherent, narrow-band excitation required for efficient and sustained VIV lock-in. This combination of intensified forcing and retained structural coherence is the physical realisation of the geometric sweet spot identified in the force data.

At Grooved AR=0.50, a qualitative breakdown in wake organisation becomes apparent. The vortex cores elongate and fragment as they convect downstream, losing the compact, well-defined structure observed in all higher aspect ratio configurations. The shear layers detach prematurely and fail to sustain fully coherent roll-up over the abbreviated streamwise body length, producing elongated, scattered vortical filaments rather than a clean alternating vortex pair. This structural fragmentation provides a direct physical explanation for the non-monotonic Cℓ,rms behaviour: at AR=0.50, the bluffness is sufficient to shorten the recirculation zone and broaden the wake, but insufficient boundary-layer development length is available for the groove perturbations to organise the shear layer into the tight, high-intensity vortex cores that characterise the AR=0.75 optimum.

3.5 Energy-Oriented Wake Assessment

3.5.1 Wake Fluctuation Energy Distribution

Fig. 14 presents the spatial distribution of the wake fluctuation energy (WFE=12⟨u′2+v′2⟩/U∞2) for the four key Stage II configurations. Here WFE denotes the kinetic energy of the resolved periodic velocity fluctuations about the time mean, used purely as a fluctuation-energy diagnostic rather than a turbulence-model quantity. The WFE field quantifies the spatial density of the aerodynamic fluctuation energy in the wake and provides a direct spatial visualisation of the geometric amplification mechanism underlying the performance hierarchy.

Figure 14: Spatial distribution of wake fluctuation energy (WFE) in the wake region Re = 200.

For Baseline AR=1.0, the WFE is distributed in elongated, parallel lobes extending far downstream along the shear-layer trajectories, with peak values concentrated approximately 4D–6D behind the cylinder base. The Grooved AR=1.0 configuration produces a slightly more compact and more intense WFE distribution compared to the smooth baseline, with the peak energy region displaced slightly upstream. This modest increase is consistent with the groove-induced acceleration of shear-layer roll-up and corresponds to the 21.1%Cℓ,rms increment identified in Section 3.2.

At Grooved AR=0.75, the WFE distribution undergoes a qualitative transformation. The high-intensity region concentrates into a wide bilobal structure immediately behind the cylinder base, with the peak WFE approximately two to three times higher than the Baseline AR=1.0 peak and located substantially closer to the body. This near-body intensification is a direct spatial consequence of the shortened recirculation zone (Lf=0.61D): coherent vortex roll-up is completed closer to the cylinder, depositing the fluctuation energy into a compact, high-density region that exerts a maximum alternating pressure force on the body. The lateral extent of the high-WFE lobes at AR=0.75 also broadens substantially relative to the circular cylinder cases, reflecting the wider separated shear layers produced by the vertically-oriented elliptical profile. Together, the near-body location, the elevated peak intensity, and the broader lateral extent of the WFE field at Grooved AR=0.75 provide the spatial energy signature of the aerodynamically optimum.

At Grooved AR=0.50, the WFE peak moves even closer to the body (Lf=0.56D), and the wake broadens further in the lateral direction. However, the spatial organisation of the WFE field becomes less structured compared to AR=0.75: the bilobal high-energy region diffuses more rapidly in the downstream direction, consistent with the vortex fragmentation and modal incoherence documented in Section 3.3. This WFE diffusion means that a fraction of the available fluctuation energy at AR=0.50 is distributed into incoherent spatial scales that cannot be efficiently captured by a narrow-band structural resonator, providing the energy-field-level explanation for why Grooved AR=0.50 underperforms Grooved AR=0.75 despite greater geometric bluffness.

3.5.2 Projected Energy Harvesting Amplification

Table 7 presents the complete aerodynamic and projected aerodynamic excitation metrics for all six Stage II configurations. The projected extractable power Pest is computed using the 1-DOF structural projection model (Eq. (8)) with representative structural parameters m∗=2.0 and ζ=0.05, which are representative of the low-mass-ratio and low-total-damping conditions typically employed in hydrokinetic VIV harvesters, where large-amplitude lock-in oscillations are commonly observed [49–51]. Because the mass, damping ratio and natural frequency are held identical across all configurations, they cancel in every inter-configuration comparison, so that the amplification factor of any configuration relative to the smooth baseline depends only on the square of the ratio of their root-mean-square lift coefficients. The selected structural parameters therefore set only the absolute magnitude of Pest and have no influence on the relative ranking of the configurations.

Fig. 15 illustrates the estimated harvested power and corresponding amplification factor for all six configurations relative to the smooth AR=1.0 baseline. The figure clearly reveals three distinct physical regimes. In the first regime, represented by the AR=1.0 pair, groove application alone raises Pest from 0.2313 to 0.3873 W/m (1.67-fold), demonstrating that the inverted biomimetic groove paradigm is effective even on the canonical circular geometry. In the second regime, the intermediate bluffness of AR=0.75 produces the optimum: Grooved AR=0.75 achieves Pest=0.5184 W/m, a 2.24-fold amplification that exceeds all other configurations. In the third regime, the maximum bluffness of AR=0.50 produces a diminishing return: despite the higher drag penalty (Cd¯=1.987), the Grooved AR=0.50 achieves only 1.89-fold amplification, decisively below the AR=0.75 result.

Figure 15: Projected aerodynamic excitation power and normalized energy amplification, referenced to the smooth circular baseline (AR=1.0), across the evaluated configurations Re = 200.

The counter-intuitive underperformance of Grooved AR=0.50 relative to Grooved AR=0.75 is the central finding of this study and is fully explained by the convergent evidence from all five analysis tools. From the force data (Fig. 7), the groove Cℓ,rms increment collapses from 13.0% at AR=0.75 to 7.0% at AR=0.50. From the spectral analysis (Fig. 9), the groove has no influence on the Strouhal number at AR=0.50, indicating that the groove perturbation cannot modify the shedding dynamics at this bluffness level. From the mean wake topology (Fig. 12), the recirculation zone at AR=0.50 is so compressed (Lf=0.56D) that insufficient boundary-layer development length remains for the groove ridges to coherently trip the shear layer. From the instantaneous vorticity (Fig. 13), the vortex roll-up at AR=0.50 is spatially fragmented and structurally incoherent. From the POD analysis (Fig. 10), the two-mode energy fraction drops from 96.2% at AR=0.75 to 88.1% at AR=0.50, and from the WFE maps (Fig. 14), the spatial energy at AR=0.50 diffuses downstream into incoherent scales unsuitable for narrow-band mechanical coupling. Together, these five lines of evidence form a physically coherent and self-consistent narrative: the Grooved AR=0.75 geometry represents the optimal balance point in a three-way competition between geometric bluffness for elevated force magnitude, groove-induced shear-layer tripping for enhanced coherence, and boundary-layer development length for sustained micro-vortex interaction.

This optimum is established among the three aspect ratios evaluated (AR=1.0, 0.75, and 0.50). The non-monotonic variation of Cℓ,rms and Pest, which rise from AR=1.0 to a maximum at AR=0.75 and then fall at AR=0.50, brackets an optimum in the neighbourhood of AR=0.75, whose precise location would require a finer aspect-ratio sweep. Accordingly, AR=0.75 is reported as the best-performing configuration among the geometries tested rather than as a globally optimal aspect ratio.

To quantify the proposed boundary-layer interaction mechanism, the key wake and boundary-layer diagnostics are consolidated in Table 8 for the four configurations: the normalised recirculation length Lf/D, the base suction −C¯p,base (the magnitude of the time-mean pressure coefficient at the rear stagnation point), the boundary-layer thickness relative to the groove amplitude δ/Ag (with Ag=D⋅AR⋅Amp and δ from the laminar boundary-layer scaling), and the two-mode POD energy fraction EPOD as a measure of shedding coherence. While Lf/D contracts and δ/Ag grows monotonically with bluffness, both Cℓ,rms and −C¯p,base peak at AR=0.75 and decline at AR=0.50, where EPOD also drops sharply from 96.2% to 88.1%.

This non-monotonic response of the aerodynamic indicators to monotonically increasing geometric bluffness suggests that the groove-organised shedding saturates beyond AR=0.75: further bluffness reduction to AR=0.50 does not strengthen it, identifying AR=0.75 as the best-performing configuration among those evaluated.

3.6 Effect of Reynolds Number on Geometric Optimum

The analysis presented in Sections 3.1–3.5 establishes Grooved AR=0.75 as the geometric optimum at Re=200. A natural question arises as to whether this optimum persists across the laminar regime or whether it is confined to the specific Reynolds number examined. To address this question, the two Grooved configurations identified as the leading candidates for energy harvesting (AR=0.75 and AR=0.50) and their corresponding smooth counterparts are re-evaluated in Stage III at two additional Reynolds numbers, Re=100 and Re=150, that span the lower half of the laminar deployment window relevant to ultra-low-speed hydrokinetic micro-harvesters. The resulting aerodynamic metrics and the corresponding groove contribution, ΔCℓ,rms, defined as the percentage increment of the grooved case over its smooth counterpart at the same aspect ratio and Reynolds number, are summarised in Table 9.

Two complementary observations emerge from Table 9. First, the groove contribution ΔCℓ,rms at AR=0.75 remains narrowly confined to the 13.0%–14.7% band across the full Reynolds range, whereas the corresponding contribution at AR=0.50 drops from 9.22% at Re=100 to 6.98% at Re=150 and 7.00% at Re=200. The saturation mechanism governing the groove effectiveness at AR=0.50 is therefore not an artefact of the Reynolds number selected in the main analysis but a robust geometric feature of the laminar regime. Second, despite this uniformly stronger groove response at AR=0.75, the absolute Cℓ,rms does not rank the two configurations identically at all Reynolds numbers. At Re=100, the absolute excitation is higher for Grooved AR=0.50 (0.4074) than for Grooved AR=0.75 (0.3610); at Re=150, the two configurations produce nearly identical values (0.5687 vs. 0.5645); and at Re=200, the ordering reverses, with Grooved AR=0.75 (0.7134) clearly outperforming Grooved AR=0.50 (0.6555). This Reynolds-dependent ordering is visualised in Fig. 16.

Figure 16: Root-mean-square lift coefficient (Cℓ,rms) as a function of Reynolds number for grooved AR=0.75 and grooved AR=0.50.

The physical mechanism underlying this crossover is clarified by the instantaneous vorticity and WFE distributions presented in Figs. 17 and 18. At Re=100 (upper rows), both configurations sustain a well-organised Kármán vortex street, with compact alternating vortex cores convecting smoothly downstream. Grooved AR=0.50 produces visibly larger and more intense vortex cores than Grooved AR=0.75, consistent with its stronger absolute bluffness and higher Cℓ,rms at this Reynolds number. The corresponding WFE distribution at Re=100 shows a broader and moderately more energetic wake for Grooved AR=0.50, fully consistent with the absolute force ranking. In this low-Reynolds regime, the boundary layer is sufficiently thick that the groove ridges interact coherently with the separating shear layer at both aspect ratios, and the unmodified bluffness advantage of AR=0.50 translates directly into a higher alternating lift magnitude.

Figure 17: Instantaneous normalised vorticity fields at Re = 100 (upper row) and Re = 150 (lower row) for grooved AR=0.75 and grooved AR=0.50. The vorticity is normalised as ωz∗=ωzD/U∞.

Figure 18: Spatial distribution of wake fluctuation energy (WFE) at Re=100 (upper row) and Re=150 (lower row) for grooved AR=0.75 and grooved AR=0.50. The corresponding Re=200 distributions are shown in Fig. 14.

A qualitatively different picture emerges at Re=150. Grooved AR=0.75 continues to exhibit the clean, alternating Kármán structure observed at Re=100, with compact vortex cores and regular downstream convection. In contrast, the Grooved AR=0.50 wake begins to lose its structural integrity: the vortex cores elongate, the alternating spacing becomes irregular, and early signs of the spatial fragmentation documented at Re=200 in Section 3.4 start to appear. The WFE field at Re=150 reflects this transition: although the two peak intensities are comparable, the AR=0.50 distribution exhibits a downstream diffusion tail that is absent from the AR=0.75 case. These qualitative changes at Re=150, together with the near-identical Cℓ,rms values of the two configurations, indicate the onset of the boundary-layer saturation mechanism at AR=0.50, partially offsetting its bluffness advantage.

Beyond this point, the saturation mechanism becomes fully developed. At Re=200, the wake at AR=0.50 is spatially fragmented, modally incoherent (cumulative two-mode energy fraction of 88.1%), and characterised by a diffused WFE field, as documented in Sections 3.3–3.5. Simultaneously, the AR=0.75 configuration retains the compact near-body WFE footprint and the high modal coherence that together define the geometric optimum. The Reynolds-dependent crossover is therefore not a coincidence but a direct consequence of the progressive thinning of the laminar boundary layer with increasing Re: once the boundary layer becomes sufficiently thin, the compressed streamwise extent of the AR=0.50 body can no longer accommodate coherent groove-to-shear-layer interaction, and the groove effectiveness collapses relative to AR=0.75.

This boundary-layer mechanism can be expressed more quantitatively. Using the classical laminar boundary-layer scaling δ99(x)=5νx/U∞ [52] and taking the streamwise body extent as x=b=AR⋅D, the normalized boundary-layer thickness is expected to vary approximately with AR/Re. Therefore, the boundary layer for the AR=0.50 body is approximately 29% thinner at Re = 200 than at Re = 100 (since 100/200=0.71). As the boundary layer thins, the compressed streamwise extent of the AR=0.50 body provides progressively less development length for coherent interaction with the groove ridges, whose amplitude relative to D is fixed at Amp⋅AR=0.025. This trend is consistent with the measured collapse of the groove contribution at AR=0.50 from 9.22% at Re = 100 to 7.0% at Re = 200, whereas the longer AR=0.75 body appears to retain a boundary layer sufficiently thick relative to the ridge scale to sustain coherent tripping, with its groove contribution remaining between 12.99% and 14.68% across the Reynolds range investigated in the present study.

The combined evidence from Table 9, Fig. 16, and the vorticity and WFE fields of Figs. 17 and 18 leads to a clear and physically coherent conclusion. The saturation of the groove contribution at AR=0.50 is a robust feature of the laminar regime and manifests at every Reynolds number examined. However, the absolute geometric optimum for VIV energy harvesting is Reynolds-number dependent: AR=0.50 delivers the stronger excitation at Re=100, the two configurations converge at Re=150, and AR=0.75 takes over as the optimum at Re=200 and, by extrapolation of the saturation trend, throughout the upper half of the laminar regime. For the hydrokinetic micro-harvester deployment conditions that motivated the present study, in which the operating Reynolds number typically lies in the upper laminar range, the Grooved AR=0.75 configuration identified in the main analysis is therefore confirmed as the physically robust best-performing configuration among those evaluated.

4  Conclusion

This study investigated the passive geometric amplification of aerodynamic excitation forces for VIV energy harvesting through a two-stage computational framework combining Taguchi L9 groove optimisation and fixed-cylinder laminar CFD. Six configurations spanning three smooth baselines and three optimally grooved elliptical cylinders at AR=1.0, 0.75, and 0.50 were evaluated at Re=200 through aerodynamic force analysis, spectral characterisation, wake topology, instantaneous vorticity, POD modal decomposition, and a 1-DOF structural projection model. The analysis was subsequently extended to Re=100 and Re=150 to establish the robustness of the identified geometric optimum across the laminar deployment window.

The Taguchi L9 analysis conclusively identified groove amplitude as the dominant control parameter, with the influence of wave number remaining negligible within the evaluated range. The optimal groove (N=24, Amp=5%) delivers a 21% increase in Cℓ,rms and a 12.5% downward Strouhal shift on the circular baseline, broadening the structural lock-in range for hydrokinetic micro-harvester deployment at ultra-low flow velocities.

Contrary to the conventional role of biomimetic grooves as drag reduction and vibration suppression devices, the present study demonstrates that when biomimetic grooves are combined with an optimally bluff elliptical geometry, they can amplify coherent aerodynamic excitation and increase the projected aerodynamic excitation power relevant to VIV energy harvesting. The most significant and physically novel finding of the Re=200 analysis is the identification of a geometric sweet spot at AR=0.75. The Grooved AR=0.75 configuration achieves the highest projected aerodynamic excitation power (Pest=0.5184 W/m), corresponding to a 2.24-fold amplification relative to the smooth circular baseline and outperforming the more aggressively bluff Grooved AR=0.50 configuration (1.89-fold). This counter-intuitive result is not a numerical artefact but a physically coherent outcome supported by five independent lines of evidence: the groove Cℓ,rms increment collapses from 13.0% to 7.0% as AR decreases from 0.75 to 0.50; groove-induced Strouhal modification disappears at AR=0.50; the recirculation zone at AR=0.50 is too compressed (Lf=0.56D) to permit coherent groove-shear-layer interaction; the instantaneous vortex roll-up at AR=0.50 fragments spatially; and the POD two-mode energy fraction drops from 96.2% (AR=0.75) to 88.1% (AR=0.50), with anomalous energy leakage into higher-order incoherent structures. Together, these observations confirm that at AR=0.50 the shortened streamwise body dimension suppresses the groove tripping mechanism, establishing AR=0.75 the best-performing interaction case among the evaluated aspect ratios at Re=200.

The Reynolds-number sensitivity analysis conducted at Re=100 and 150 refines and generalises this finding. The groove contribution at AR=0.75 remains narrowly confined to the 12.99%–14.68% band across the full Reynolds range, whereas the corresponding contribution at AR=0.50 drops from 9.22% at Re=100 to approximately 7% at Re=150 and Re=200. The boundary-layer saturation mechanism that suppresses the groove effectiveness at AR=0.50 is therefore a robust geometric feature of the laminar regime rather than a Reynolds-specific artefact. At the same time, the absolute geometric optimum is Reynolds-number dependent: AR=0.50 delivers the stronger excitation at Re=100 on account of its unmodified bluffness advantage in the thicker boundary-layer regime, the two configurations converge at Re=150, and AR=0.75 emerges as the optimum at Re=200 once the boundary layer has thinned sufficiently for the saturation mechanism to be fully activated. This indicates that the geometric optimum shifts from AR=0.50 to AR=0.75 at Re=150, corresponding to the onset of the boundary-layer saturation mechanism. For hydrokinetic micro-harvester deployment conditions in the upper laminar range, the Grooved AR=0.75 configuration is therefore confirmed as the physically robust geometric optimum.

A further structural advantage of the Grooved AR=0.75 configuration is the restoration of the Strouhal number from St=0.2197 (smooth AR=0.75) to St=0.1953 following groove application, aligning the shedding frequency with the circular cylinder reference value. A VIV harvester resonator tuned for the standard circular geometry can therefore be deployed directly with the Grooved AR=0.75 body without structural retuning, a practical advantage for scalable device implementation.

The present study is subject to several limitations. The simplified 1-DOF projection model uses constant mass and damping parameters and does not resolve the coupled fluid-structure dynamics. Future work should therefore extend the framework to fully coupled simulations that quantify lock-in, added mass, and reduced-velocity effects on the actual energy conversion efficiency. A finer aspect-ratio sweep near AR=0.75 and a co-optimisation of the groove parameters at each aspect ratio and Reynolds number would refine the identified geometric optimum. Three-dimensional simulations would clarify the geometric synergy and whether geometry-dependent differences in the transition threshold could shift the absolute optimum. A finer Reynolds-number sweep between Re = 150 and Re = 200 would localise the crossover Reynolds number. Experimental validation in a water-channel facility would consolidate the practical viability of the proposed concept.

Acknowledgement: Not applicable.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Yunus Celik; methodology, Yunus Celik and Burhan Necati Kiziloglu; software, Yunus Celik; validation, Yunus Celik and Burhan Necati Kiziloglu; formal analysis, Yunus Celik; investigation, Yunus Celik; resources, Yunus Celik; data curation, Yunus Celik; writing—original draft preparation, Yunus Celik; writing—review and editing, Yunus Celik and Burhan Necati Kiziloglu; visualization, Yunus Celik; supervision, Yunus Celik; project administration, Yunus Celik. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available from the Corresponding Author, Yunus Celik, upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

Nomenclature

The following symbols and abbreviations are used in this manuscript:

References

1. Bernitsas MM, Raghavan K, Ben-Simon Y, Garcia EMH. VIVACE (vortex induced vibration aquatic clean energya new concept in generation of clean and renewable energy from fluid flow. J Offshore Mech Arct Eng. 2008;130(4):041101. doi:10.1115/1.2957913. [Google Scholar] [CrossRef]

2. Zhao M. A review of recent studies on the control of vortex-induced vibration of circular cylinders. Ocean Eng. 2023;285:115389. doi:10.1016/j.oceaneng.2023.115389. [Google Scholar] [CrossRef]

3. Bakhtiar S, Khan FU, Fu H, Hajjaj AZ, Theodossiades S. Fluid flow-based vibration energy harvesters: a critical review of state-of-the-art technologies. Appl Sci. 2024;14(23):11452. doi:10.3390/app142311452. [Google Scholar] [CrossRef]

4. Cao D, He J, Zeng H, Zhu Y, Chan SZ, Williams MR, et al. A review of oscillators in hydrokinetic energy harnessing through vortex-induced vibrations. Fluids. 2025;10(4):78. doi:10.3390/fluids10040078. [Google Scholar] [CrossRef]

5. Williamson CHK, Govardhan R. Vortex-induced vibrations. Annu Rev Fluid Mech. 2004;36(1):413–55. doi:10.1146/annurev.fluid.36.050802.122128. [Google Scholar] [CrossRef]

6. Sarpkaya T. A critical review of the intrinsic nature of vortex-induced vibrations. J Fluids Struct. 2004;19(4):389–447. doi:10.1016/j.jfluidstructs.2004.02.005. [Google Scholar] [CrossRef]

7. Khalak A, Williamson CHK. Motions, forces and mode transitions in vortex-induced vibrations at low mass-damping. J Fluids Struct. 1999;13(7–8):813–51. doi:10.1006/jfls.1999.0236. [Google Scholar] [CrossRef]

8. Wang J, Geng L, Ding L, Zhu H, Yurchenko D. The state-of-the-art review on energy harvesting from flow-induced vibrations. Appl Energy. 2020;267(15):114902. doi:10.1016/j.apenergy.2020.114902. [Google Scholar] [CrossRef]

9. Moradi MA, Mojra A. Flow and noise control of a cylinder using grooves filled with porous material. Phys Fluids. 2024;36(4):045134. doi:10.1063/5.0205125. [Google Scholar] [CrossRef]

10. Zhao F, Zeng L, Bai H, Alam MM, Wang Z, Dong Y, et al. Vortex-induced vibration of a sinusoidal wavy cylinder: the effect of wavelength. Phys Fluids. 2024;36(8):081909. doi:10.1063/5.0219753. [Google Scholar] [CrossRef]

11. Talley S, Mungal G, Iaccarino G. Flow around cactus-shaped cylinders. In: Annual research briefs. Stanford, CA, USA: Center for Turbulence Research, Stanford University; 2001. p. 363–74. [Google Scholar]

12. Zdravkovich MM. Flow around circular cylinders. Vol. 1, Fundamentals. New York, NY, USA: Oxford University Press; 1997. [Google Scholar]

13. Lam K, Lin YF. Effects of wavelength and amplitude of a wavy cylinder in cross-flow at low Reynolds numbers. J Fluid Mech. 2009;620:195–220. doi:10.1017/S0022112008004217. [Google Scholar] [CrossRef]

14. Babu P, Mahesh K. Aerodynamic loads on cactus-shaped cylinders at low Reynolds numbers. Phys Fluids. 2008;20(3):035112. doi:10.1063/1.2887982. [Google Scholar] [CrossRef]

15. Talley S, Iaccarino G, Mungal G, Mansour NN. An experimental and computational investigation of flow past cacti. In: Annual research briefs. Stanford, CA, USA: Center for Turbulence Research, NASA Ames/Stanford University; 2001. p. 51–63. [Google Scholar]

16. Williamson CHK. Vortex dynamics in the cylinder wake. Annu Rev Fluid Mech. 1996;28(1):477–539. doi:10.1146/annurev.fl.28.010196.002401. [Google Scholar] [CrossRef]

17. Shi X, Alam MM, Bai H. Wakes of elliptical cylinders at low Reynolds number. Int J Heat Fluid Flow. 2020;82(1):108553. doi:10.1016/j.ijheatfluidflow.2020.108553. [Google Scholar] [CrossRef]

18. Park J, Kwon K, Choi H. Numerical solutions of flow past a circular cylinder at Reynolds numbers up to 160. KSME Int J. 1998;12(6):1200–5. doi:10.1007/BF02942594. [Google Scholar] [CrossRef]

19. Jackson CP. A finite-element study of the onset of vortex shedding in flow past variously shaped bodies. J Fluid Mech. 1987;182:23–45. doi:10.1017/S0022112087002234. [Google Scholar] [CrossRef]

20. Shi X, Alam MM, Zhu H, Ji C, Bai H, Sharifpur M. Flow three-dimensionality of wavy elliptic cylinder: vortex shedding bifurcation. Ocean Eng. 2024;301(1):117527. doi:10.1016/j.oceaneng.2024.117527. [Google Scholar] [CrossRef]

21. Pradhan A, Arif MR, Afzal MS, Gazi AH. On the origin of forces in the wake of an elliptical cylinder at low Reynolds number. Environ Fluid Mech. 2022;22(6):1307–31. doi:10.1007/s10652-022-09892-z. [Google Scholar] [CrossRef]

22. Kim J, Choi H. Distributed forcing of flow over a circular cylinder. Phys Fluids. 2005;17(3):033103. doi:10.1063/1.1850151. [Google Scholar] [CrossRef]

23. Bearman PW. Vortex shedding from oscillating bluff bodies. Annu Rev Fluid Mech. 1984;16(1):195–222. doi:10.1146/annurev.fl.16.010184.001211. [Google Scholar] [CrossRef]

24. Law YZ, Jaiman RK. Passive control of vortex-induced vibration by spanwise grooves. J Fluids Struct. 2018;78(4):297–315. doi:10.1016/j.jfluidstructs.2018.08.004. [Google Scholar] [CrossRef]

25. Huang S. VIV suppression of a two-degree-of-freedom circular cylinder and drag reduction of a fixed circular cylinder by the use of helical grooves. J Fluids Struct. 2011;27(7):1124–33. doi:10.1016/j.jfluidstructs.2011.07.005. [Google Scholar] [CrossRef]

26. Zhou B, Wang X, Guo W, Gho WM, Tan SK. Experimental study on flow past a circular cylinder with rough surface. Ocean Eng. 2015;109(Pt. 1):7–13. doi:10.1016/j.oceaneng.2015.08.062. [Google Scholar] [CrossRef]

27. Han P, Huang Q, Pan G, Qin D, Wang W, Goncalves RT, et al. Optimal energy harvesting efficiency from vortex-induced vibration of a circular cylinder. Ocean Eng. 2023;282(2):114869. doi:10.1016/j.oceaneng.2023.114869. [Google Scholar] [CrossRef]

28. Fershalov A, Elvin N, Orlandini P, Avros I, Liu Y. Harvesting sustainable energy through vortex-induced vibrations of finite length cylinder. Phys Fluids. 2025;37(4):044108. doi:10.1063/5.0260317. [Google Scholar] [CrossRef]

29. Zhang M, Wu T, Xu F. Vortex-induced vibration of bridge decks: describing function-based model. J Wind Eng Ind Aerodyn. 2019;195(1368):104016. doi:10.1016/j.jweia.2019.104016. [Google Scholar] [CrossRef]

30. Morse TL, Williamson CHK. Prediction of vortex-induced vibration response by employing controlled motion. J Fluid Mech. 2009;634:5–39. doi:10.1017/S0022112009990516. [Google Scholar] [CrossRef]

31. Zhang M, Zhang C, Abdelkefi A, Yu H, Gaidai O, Qin X, et al. Piezoelectric energy harvesting from vortex-induced vibration of a circular cylinder: effect of Reynolds number. Ocean Eng. 2021;235(1):109378. doi:10.1016/j.oceaneng.2021.109378. [Google Scholar] [CrossRef]

32. Barkley D, Henderson RD. Three-dimensional Floquet stability analysis of the wake of a circular cylinder. J Fluid Mech. 1996;322:215–41. doi:10.1017/S0022112096002777. [Google Scholar] [CrossRef]

33. Hasheminejad SM, Jarrahi M. Numerical simulation of two dimensional vortex-induced vibrations of an elliptic cylinder at low Reynolds numbers. Comput Fluids. 2015;107(1):25–42. doi:10.1016/j.compfluid.2014.10.011. [Google Scholar] [CrossRef]

34. Hu Y, Yuan H, Shu S, Niu X, Li M. An improved momentum exchanged-based immersed boundary–lattice Boltzmann method by using an iterative technique. Comput Math Appl. 2014;68(3):140–55. doi:10.1016/j.camwa.2014.05.013. [Google Scholar] [CrossRef]

35. Rajani BN, Kandasamy A, Majumdar S. Numerical simulation of laminar flow past a circular cylinder. Appl Math Model. 2009;33(3):1228–47. doi:10.1016/j.apm.2008.01.017. [Google Scholar] [CrossRef]

36. Celik Y. Surrogate-based robust design of asymmetric splitter plates for aerodynamic performance efficiency of flatback airfoils. Phys Fluids. 2026;38(4):047108. doi:10.1063/5.0325553. [Google Scholar] [CrossRef]

37. Celik Y, Ingham D, Ma L, Kiziloglu BN, Pourkashanian M. Numerical analysis of unsteady performance of an aerofoil with surface openings under Darrieus motion. Wind Struct. 2026;42(3):329–55. doi:10.12989/was.2026.42.3.329. [Google Scholar] [CrossRef]

38. Taira K, Brunton SL, Dawson STM, Rowley CW, Colonius T, McKeon BJ, et al. Modal analysis of fluid flows: an overview. AIAA J. 2017;55(12):4013–41. doi:10.2514/1.J056060. [Google Scholar] [CrossRef]

39. Posdziech O, Grundmann R. A systematic approach to the numerical calculation of fundamental quantities of the two-dimensional flow over a circular cylinder. J Fluids Struct. 2007;23(3):479–99. doi:10.1016/j.jfluidstructs.2006.09.004. [Google Scholar] [CrossRef]

40. Henderson RD. Nonlinear dynamics and pattern formation in turbulent wake transition. J Fluid Mech. 1997;352:65–112. doi:10.1017/S0022112097007465. [Google Scholar] [CrossRef]

41. Williamson CHK. Oblique and parallel modes of vortex shedding in the wake of a circular cylinder at low Reynolds numbers. J Fluid Mech. 1989;206:579–627. doi:10.1017/S0022112089002429. [Google Scholar] [CrossRef]

42. Taguchi G. System of experimental design: engineering methods to optimize quality and minimize costs. Vol. 1 & 2. White plains, NY, USA: UNIPUB/Kraus International Publications; 1987. [Google Scholar]

43. Erkan O, Ozkan M, Celik Y, Khalid MSU. Taguchi-based multi-factor analysis of self-starting behavior in vertical axis wind turbine farms. Energy. 2025;341(1–3):139495. doi:10.1016/j.energy.2025.139495. [Google Scholar] [CrossRef]

44. Lumley JL. The structure of inhomogeneous turbulent flows. In: Yaglom AM, Tatarski VI, editors. Atmospheric turbulence and radio wave propagation. Moscow, Russia: Nauka; 1967. p. 166–78. [Google Scholar]

45. Berkooz G, Holmes P, Lumley JL. The proper orthogonal decomposition in the analysis of turbulent flows. Annu Rev Fluid Mech. 1993;25(1):539–75. doi:10.1146/annurev.fl.25.010193.002543. [Google Scholar] [CrossRef]

46. Sirovich L. Turbulence and the dynamics of coherent structures. Part I: coherent structures. Q Appl Math. 1987;45(3):561–71. doi:10.1090/qam/910462. [Google Scholar] [CrossRef]

47. Jie H, Liu YZ. Large eddy simulation of turbulent flow over a cactus-analogue grooved cylinder. J Vis. 2015;18(2):171–88. doi:10.1007/s12650-015-0294-x. [Google Scholar] [CrossRef]

48. Lin JC, Hsieh SC. Convection velocity of vortex structures in the near wake of a circular cylinder. J Eng Mech. 2003;129(10):1108–18. doi:10.1061/(ASCE)0733-9399(2003)129:10(1108). [Google Scholar] [CrossRef]

49. Sun H, Kim ES, Nowakowski G, Mauer E, Bernitsas MM. Effect of mass-ratio, damping, and stiffness on optimal hydrokinetic energy conversion of a single, rough cylinder in flow induced motions. Renew Energy. 2016;99(6):936–59. doi:10.1016/j.renene.2016.07.024. [Google Scholar] [CrossRef]

50. Pal A, Soti AK. Power harvesting from VIV of rigidly-coupled cylinders in tandem arrangement. In: Fluid mechanics and fluid power, Volume 5. Lecture notes in mechanical engineering. Singapore: Springer; 2023. p. 621–30. doi:10.1007/978-981-19-7055-9_53. [Google Scholar] [CrossRef]

51. Barrero-Gil A, Pindado S, Avila S. Extracting energy from vortex-induced vibrations: a parametric study. Appl Math Model. 2012;36(7):3153–60. doi:10.1016/j.apm.2011.09.085. [Google Scholar] [CrossRef]

52. Schlichting H, Gersten K. Boundary-layer theory. 9th ed. Berlin, Germany: Springer; 2017. doi:10.1007/978-3-662-52919-5. [Google Scholar] [CrossRef]