Numerical Investigation of Aeroacoustic Response of an Oscillating-Wing Power Extractor

1 Introduction

The demand for green energy is continuously increasing due to its significant role in promoting economic growth, mitigating global warming, and reducing environmental pollution. Among the various renewable energy sources, wind energy has attracted considerable interest because of its advantages, including zero carbon emissions, operational safety, and environmental sustainability [1]. Traditionally, wind energy has been harnessed using rotary turbines. Although these systems are mature and highly efficient, they are often noisy and their rotating blades are susceptible to high centrifugal and cyclic loading. Moreover, rotary wind turbines are optimized for specific operating conditions, and their aerodynamic performance drops when operating away from the optimal tip-speed ratio [2,3,4]. This has motivated the development of alternative wind energy harvesting technologies, including airborne wind turbines, floating offshore wind turbines, Savonius and Darrieus vertical-axis wind turbines, and oscillating-wing power extractors (OWPEs) [5,6].

In recent times, OWPE has emerged as a unique and promising technology that harnesses energy through the oscillatory motion of a wing, inspired by flapping motion observed in natural flyers and swimmers [7]. An example of this technology is the dual-wing generator designed, developed and tested by Festo Corporate [8]. Oscillating foil-based turbine was first suggested by Mckinney and DeLaurier [9]. Unlike rotary turbine, an oscillating foil can operate at all incoming velocities as the reduced frequency f* and pitch amplitude θo can be changed dynamically. Moreover, it is an ideal choice for small and portable applications, a looming market in the energy sector [10]. The concept of the flapping wing energy harvester has been discussed in detail in the review articles by Alam et al. [11] and Liu et al. [12]. Kinsey and Dumas [13] investigated the performance of a single oscillating turbine for 0° < θo < 90° and 0 < f* < 0.25 using 2-D numerical simulations at Re = 1100 and Ho/c = 1, where Ho and c are the heave amplitude and chord length, respectively. They have shown that efficiency of a single turbine can reach as high as 35% when θo = 75° and f* = 0.15. Further, they have categorized the entire regime into two, propulsion and power extraction using a parameter known as the feathering limit.

Previous researchers have investigated OWPE concept using experiments and numerical study. Zhu [14] investigated the wake stability and energy harvesting performance of a single oscillating wing at Re = 1000. It was shown that the power extraction is closely related with the evolution of the wake. An efficient wake evolution can lead to more power extraction and efficiency. Kinsey et al. [15] built and tested a new prototype of 2 kW oscillating foil using experiments at Re = 4.8 × 105. They have concluded that the efficiency of a single oscillating foil can reach ~35% after taking into account the mechanical losses. Ma et al. [16] investigated the energy harvesting behavior of a flapping wing in arc trajectory at Re = 1 × 106. They have shown that both arc radius R and heave amplitude Ho play a vital role on the energy harvesting characteristics of flapping wing. Recently, Sitorus and Ko [8] compared the energy harvesting characteristics of three types of flapping wing: pitch-heave, left swing and right swing at Re = 1.7 × 106. They have shown that the flapping wing in left swing outperforms pitch-heave and right swing energy harvesting characteristics.

A single oscillating foil extracts energy from the flow through its oscillatory motion. To enhance energy extraction, researchers have investigated systems comprising multiple oscillating foils, where interactions between the foils can influence the overall energy harvesting performance [17]. Kinsey and Dumas [18] investigated the performance of tandem oscillating foils using two-dimensional simulations at Ho/c = 1 and Re = 9 × 105. They demonstrated that the efficiency of such a system could reach as high as 64%. In tandem configurations, interactions between the foils can be either favorable or unfavorable, resulting in an increase or decrease in the overall energy harvesting efficiency, respectively. Karbasian et al. [19] investigated the power extraction performance of multiple flapping foils in tandem formation at Ho/c = 1 and Re = 9 × 105. They found that the power extraction performance drops after the second oscillating foil. Xu et al. [20] measured the power extraction performance of a tandem oscillating foil using experiments at Re = 2.6 × 105–3.8 × 105. They showed that the optimal reduced frequency for achieving highest efficiency of 29% is 0.1. Ashraf et al. [21] numerically studied the performance of single and tandem oscillating wings at Ho/c = 1 and Re = 2 × 104. Their results showed that the efficiency of a single wing could reach approximately 34% for f* = 0.13 and θo = 73°. However, the performance of the tandem wing configuration was observed to be around 54% when the wings oscillated in opposite directions and were separated by a distance of 2c.

A few other bio-inspired concepts and power enhancement techniques have also been proposed to enhance the overall performance of oscillating-wing energy harvesters. Le et al. [22] elucidated the characteristics of energy harvester inspired from Scallops-shell. They found that the size and strength of LEV are strongly affected by morphology of the mimicked wings, the effects of which can lead to improvement in the energy harvesting efficiency. Lahooti and Kim [23] investigated the influence of an upstream body on the power extraction performance of an oscillating wing at Re = 1000. By placing the upstream body at a proper location, the power extraction efficiency of an oscillating wing can be increased by 30%.

One of the major challenges of wind turbines is the aerodynamic noise generated by these systems, particularly when they are deployed near residential and urban environments. Several studies have investigated the aerodynamic noise generated by conventional wind turbines using numerical simulations [24,25,26]. A few other studies have evaluated the aero-acoustics noise of Darrieus wind turbines [27,28]. These studies have shown that aerodynamic noise is strongly influenced by flow structures generated around the energy extraction system. Unlike conventional wind turbines, OWPEs extract energy through periodic pitching and heaving motions of a wing. Due to their compact size, OWPEs are suitable for installation on residential buildings, rooftops, and other urban environments. However, the unsteady aerodynamic mechanisms associated with oscillating wings generate complex vortex structures that may significantly influence aeroacoustic characteristics. Consequently, understanding the aerodynamic noise generated by OWPEs is essential for assessing their practical feasibility in urban environments.

Against this backdrop, the present study investigates the aerodynamic noise characteristics of an OWPE and examines the effects of key wing kinematic parameters on noise generation. Although the aerodynamic performance and energy extraction characteristics of OWPEs have been investigated in several previous studies, their aeroacoustic behavior has not yet been investigated. To the best of the authors’ knowledge, no previous study has systematically examined the aerodynamic noise generated by OWPEs or evaluated the influence of wing kinematics on aerodynamic noise generation. This lack of knowledge represents an important research gap, especially because OWPEs are expected to be installed in locations where noise is an important consideration. In addition, a TOPSIS-based multi-criteria analysis is performed to identify the optimal operating condition by simultaneously considering the power extraction performance and acoustic levels. The findings of this work can provide some insights into the aeroacoustic behavior of oscillating-wing energy harvesters. Furthermore, it can contribute to the development of quieter OWPE systems suitable for deployment in residential and urban environments.

2 Methodology

2.1 Wing Kinematics

Table 1 presents the geometric and operating parameters of OWPE used in this study. The airfoil was modeled as a NACA0015 airfoil. These parameters define the wing motion and the conditions considered for the analysis.

The wing kinematics of the OWPE can be represented mathematically by simultaneous heave motion H(t) and pitch motion θ(t) as given by Eqs. (1) and (2). The pitch center XP is located at c/3 from the leading edge (LE). Ht=Hosin(ωt+φ)(1) θt=θosin(ωt)(2) here, Ho/c denotes the heave amplitude (=1), f the oscillation frequency, θo the pitch amplitude, φ the heave-pitch phase shift (φ = 90°), and ω = 2πf. The Reynolds number Re and reduced frequency f* can be computed from Eqs. (3) and (4). Re=ρU∞cμ(3) f*=fcU∞(4) here, ρ denotes the fluid density, U∞ the free-stream velocity, c the chord length and μ the dynamic viscosity of fluid. The aerodynamic coefficients per unit span (CX, CY and CM) can be computed from Eqs. (5)–(7). CX=FX0.5ρU∞2c(5) CY=FY0.5ρU∞2c(6) CM=M0.5ρU∞2c2(7) here, FX and FY denote the horizontal and vertical forces per unit span acting on the oscillating foil, respectively, and M is the pitching moment per unit span. The instantaneous and cycle-averaged power coefficients per unit span are calculated using Eqs. (8)–(10). Pt=FYtH•t+Mtθ•t=FYtVYt+MtΩt(8) CP=P0.5ρU∞3c(9) CP¯=CP,Y¯+CP,θ¯=∫0TCYtVYtdtTU∞+∫0TCMtcΩtdtTU∞(10) here, P the instantaneous power extracted from the oscillating foil per unit span, CP the power coefficient, CP,Y the heave power coefficient, CP,θ the pitch power coefficient, V Y = H • is the heaving velocity, and Ω = θ • is the pitching angular velocity.

2.2 Numerical Methodology

Simulations are performed using the finite volume–based commercial CFD software ANSYS Fluent. The fluid is assumed to be incompressible, Newtonian, and isothermal with constant physical properties. The flow is modeled as transient due to the periodic flapping motion of the wing. The gravitational and thermal effects are neglected. The wing undergoes a prescribed sinusoidal motion with fixed amplitude and frequency, independent of the aerodynamic loads. A pressure-based transient solver is employed to solve the two-dimensional incompressible unsteady Reynolds-Averaged Navier–Stokes (URANS) equations, defined by Eqs. (11)–(13). To accommodate the mesh deformation caused by the oscillatory motion of the wing, the local remeshing option available in the dynamic mesh framework is employed to continuously adjust the mesh cells in the external flow-field region. For the numerical solution procedure, the pressure-implicit with splitting of operators (PISO) algorithm is used for pressure–velocity coupling. The PISO algorithm is selected because of its improved pressure–velocity coupling accuracy and numerical stability for transient unsteady flow simulations involving moving boundaries. A second-order spatial discretization scheme is adopted for pressure, while the momentum and turbulence equations are discretized using the second-order upwind scheme to ensure improved solution accuracy. Furthermore, a first-order implicit scheme is employed for the temporal discretization of all governing equations to ensure numerical stability during transient simulations. The convergence criterion for the continuity, momentum, and turbulence equations is set to a residual value of 10−6 at each time step.

Continuity equation: ∂ui∂xi=0(11) where ui is the mean velocity component.

Momentum equation: ∂∂tρui+∂∂xjρuiuj=−∂p∂xi+∂∂xjμ∂ui∂xj+∂uj∂xi−23δij∂ul∂xl+∂∂xj−ρui'uj'¯(12) −ρui'uj'¯=μt∂ui∂xj+∂uj∂xi−23ρk+μt∂uk∂xkδij(13) where ρ is the fluid density, uj and ul are the mean velocity components, p is the pressure, μ is the molecular viscosity, μt is the turbulent viscosity, δij is the Kronecker delta, k is the turbulent kinetic energy, and − ρ u i ′ u j ′ ¯ denotes the Reynolds stress tensor. The turbulent viscosity (μt) for the k–ω shear stress transport (SST) turbulence model, given by Eq. (14), is defined as μt=ρkω1max1α*,SF2a1ω(14) where S is the strain rate magnitude, ω is the specific dissipation rate, F2 is the blending function, and a1 and α* are the SST model constants/functions.

k–ω SST turbulence model: ∂ρk∂t+∂∂xiρkui=∂∂xjΓk∂k∂xj+Gk−Yk+Gb(15) ∂ρω∂t+∂∂xi=∂∂xjΓω∂ω∂xj+Gω−Yω+Dω+Gωb(16) where Gk and Gω are the generation terms of k and ω, Γk and Γω are the corresponding effective diffusivity, Yk and Yω are the dissipation terms for k and ω due to turbulence, Gb and Gωb account for buoyancy terms of k and ω, and Dω is the cross-diffusion term, respectively. The effective diffusion coefficients Γk and Γω are expressed as Γk=μ+μtσk(17) Γω=μ+μtσω(18) where σk and σω are the turbulent Prandtl numbers for k and ω, respectively.

The k–ω SST turbulence model, defined by Eqs. (14)–(18), is employed for turbulence closure. High numerical accuracy is essential in the present study because the unsteady pressure fluctuations developed on the wing surface are directly used as input for the FW–H acoustic model to predict the aerodynamic noise generated by the oscillating wing. Therefore, accurate prediction of the unsteady flow field around the wing is important for obtaining reliable aeroacoustic results. The k–ω SST turbulence model is adopted because it provides improved prediction of adverse pressure-gradient flows and flow separation. The model combines the near-wall accuracy of the Wilcox k–ω formulation with the robustness of the k–ε model in the outer flow region through a blending-function approach [29]. Consequently, it is capable of accurately predicting the onset and development of flow separation under unsteady flow conditions. This turbulence model is particularly suitable for flapping wing energy harvesting applications, where the periodic motion of the wing generates complex unsteady flow structures such as leading-edge vortices (LEVs), trailing-edge vortices (TEVs), and wake regions. These vortical structures strongly influence the instantaneous aerodynamic forces, power extraction performance, and aerodynamic noise characteristics of the system [30]. Furthermore, the interaction between the oscillating wing can cause significant boundary-layer modification, local flow separation, and transient wake development. Therefore, accurate modeling of near-wall turbulence and transient flow behavior is necessary to capture the aerodynamic and aeroacoustic characteristics of the oscillating-wing power extractor more reliably [31,32,33].

2.3 Computational Domain and Boundary Conditions

The details of computational domain, mesh model, and boundary conditions used in this work are illustrated in Fig. 1. The X and Y axis are aligned with the streamwise and cross-streamwise directions. The origin is fixed at the center of the wing. The computational domain is a rectangle of length l (–20c ≤ l ≤ 40c) and breadth b (–30c ≤ b ≤ 30c). The domain is subdivided into two zones, wing zone (rigid motion) and buffer zone (deforming). These two zones are separated by a non-conformal sliding interface. The NACA 0015 wing is embedded inside the wing zone (circular shape). The wing zone and the interface executes ‘simultaneous heave and pitch motion’ while the buffer zone deforms due to the flapping movement of wing. Four distinct boundaries (inlet, outlet, side boundaries and wing) are present in the domain. The boundary conditions used in the simulations: velocity inlet u = U∞ for the left boundary, pressure outlet p = Patm (101,325 Pa) for the right boundary, and symmetry for the side boundaries. The side boundaries are located at a distance of 30c from the wing and therefore, they have negligible influence on the wake structures of oscillating wing. A no-slip wall condition u = 0, v = 0 is imposed on the wing. A medium turbulent intensity of 3% is set at the inlet and outlet boundaries while the turbulent viscosity ratio is fixed to be 0.01. The simulations are performed at a Reynolds number of Re = 8.58 × 104, based on the wing chord length and free-stream velocity, corresponding to a Mach number of Ma = 0.015. Since the Mach number is less than 0.3, the compressibility effects can be neglected.

As shown in Fig. 1, unstructured triangular cells are employed in most region of computational domain. However, structured quadrilateral cells of O-type are employed around the wing to predict the boundary layer effects more accurately. A total of 8 × 104 cells is used for the entire domain and the optimum time-step size is found to be T/1000. For the k–ω shear stress transport (SST) turbulence model to function accurately, it is necessary to satisfy the criteria of y+ = 1. Therefore, the first layer of the boundary mesh is kept at a distance of 2.3 × 10–4c from the wing.

Figure 1: (a) Wing kinematics (red-mid downstroke and blue-mid upstroke). (b) Computational domain and boundary conditions adopted in the present study.

The motion of flapping wing is provided as input via user-defined function (UDF) and the mesh motion is performed using dynamic mesh method. Spring-based smoothing and remeshing strategies are used to handle the dynamic motion of mesh. A spring constant of 0.01 is used in the smoothing method. The minimum length scale is selected from the value of minimum cell edge length, while the maximum scale length is maintained at 10 times the minimum length scale. The maximum skewness of the deformed mesh during remeshing is restricted to 0.6. In the present work, spring based technique of dynamic mesh method is employed. In the dynamic mesh method of ANSYS Fluent [34], the cells inside the domain change dynamically to accommodate the flapping motion of wing. For every time-step, the cells in the domain are updated by the solver through adding, removing or stretching. The governing equation of a general scalar of ϕ, on an control volume V, for moving boundaries in dynamic mesh in integral form can be computed from Eq. (19). ddt∫VρϕdV+∫∂Vρϕu→−ug→dA→=∫∂VΓ∇ϕdA→+∫VSϕdV(19) here ρ is the fluid density, Г is the diffusion coefficient, u → is the flow velocity vector, Sϕ is the source term of ϕ, u g → is the mesh velocity of the dynamic mesh, ∂V the boundary of the control volume. The time derivative term in Eq. (19) is discretized via second order backward scheme as follows:

There are two strategies available to tackle the cell update, smoothing and remeshing. The smoothing method operates by moving the interior and boundary nodes that absorb the movement of the rigid/deforming boundary without changing the total nodes and mesh connection. Two types of smoothing strategies are implemented in ANSYS Fluent for interior surface namely spring-based smoothing and diffusion-based smoothing. The computational effort of diffusion-based smoothing is very high in relative to spring-based smoothing but it aids to improve the mesh quality slightly [34]. Spring based smoothing is an efficient and quick way smooth any cell or face that belong to moving/deforming boundary. In this strategy, the edges between the two nodes are modeled in the form of network of springs and the movement of boundaries generate spring force. From Hooke’s Law, this spring force on a node, are given by Eqs. (21) and (22). Fi→=∑jnikijΔxj→−Δxi→(21) kij=kfacxj→−xi→(22) here Δ x i → , Δ x j → denotes the displacements of node i and its neighbor j, ni the number of nodes connected to I, kij the spring constant, and kfac is the spring constant factor. The range of spring constant factor is 0 < kfac < 1. A lower value of kfac will give a less stiff grid. A higher value of kfac will primarily deform the nodes in the vicinity of the moving body [34]. At equilibrium condition, the summation of force generated from the springs interlinked to the nodes is equal to zero. This condition provides an iterative equation as shown in Eq. (23), the solution of which can be obtained from iterative process by specifying the convergence tolerance and number of iterations. The iterative solutions stops when any of these conditions are met [34].

Remeshing is a dynamic mesh strategy available in ANSYS Fluent for handling moving boundaries. It is particularly effective for simulations involving relatively large boundary displacements. In this approach, the software automatically reconstructs the mesh whenever cells or faces exceed predefined quality criteria, such as skewness or size limits, thereby preserving mesh quality and improving solution accuracy throughout the simulation. In general, both remeshing and smoothing are adopted simultaneously as it allows large time step to generate better quality grid. Local cell remeshing option marks the interior cells that exceeds the value of user specified criteria (skewness and size in terms of length scale) for every time step and locally remeshes these cells. The grid is updated if the new cells improve the skewness and the solution is interpolated from the previous cells [34]. It is important to remember that the marking of cells for remeshing based on specified criteria (skewness and size) happens before the motion of the boundary. During the remeshing process, there is no direct control over the skewness and size criteria of the remeshing cells. Therefore, one must cautious in choosing the tight criteria for both skewness and size (length scale) to effectively control the remeshing cells.

2.4 FW-H Acoustic Model

In this work, the Ffowcs Williams−Hawkings (FW–H) acoustic analogy model [35] is employed to predict the aerodynamic noise generated by the oscillating wing. The Ffowcs Williams−Hawkings (FW–H) formulation is based on Lighthill’s acoustic analogy and is widely used to predict aerodynamic noise in the mid- and far-field. It calculates the acoustic pressure at receiver locations using the unsteady flow data obtained from the transient CFD simulation. Since the aerodynamic flow and acoustic calculations are performed separately, the FW–H method requires less computational effort than fully coupled aeroacoustic simulations while still providing accurate far-field noise predictions [35].

In the present study, the x-distance of the receivers is defined as 1c, 3c, 5c, and 10c measured from the pitch axis, while the y-distance is defined as 0c, 3c, −3c, 5c, and −5c measured from the bottom-most position of the wing, where c denotes the airfoil chord length. The wing surface is used as the acoustic source surface, and the dominant aerodynamic noise is expected to arise from loading noise associated with unsteady aerodynamic forces and periodic vortex shedding. Although the present study employs a two-dimensional URANS-based FW–H model, the actual flow around the oscillating wing is inherently three-dimensional. Therefore, three-dimensional flow structures, such as spanwise vortices, are not captured. The two-dimensional assumption may overpredict the SPL and reduce the accuracy of high-frequency noise prediction. Nevertheless, the present approach is suitable for comparing the relative aeroacoustic performance under different operating conditions. The governing equation of the FW–H model is an inhomogeneous wave equation derived from the continuity and Navier–Stokes equations [34,35]. Since the present flapping wing operates under low-Mach-number conditions, the quadrupole source term was neglected, and only the monopole and dipole source terms were considered. The FW–H equation can be expressed as follows: 1a02∂2p′∂t2−∇2p′=−∂∂xiPijnj+ρuiun−vnδf+∂∂tρ0vn+ρun−vnδf(24) Pij=pδij−μ∂ui∂xj+∂uj∂xi−23∂uk∂xkδij(25) where p′ is the acoustic pressure fluctuation, a0 is the speed of sound, Pij is the compressive stress tensor, δ(f) is the Dirac delta function, nj is the outward unit normal vector, un is the normal component of fluid velocity, vn is the normal component of surface velocity, and ρ0 is the ambient density.

In ANSYS Fluent, the FW–H acoustic analysis is performed in two steps. First, the transient CFD simulation provides the time histories of pressure, velocity, and density on the wing surface, which is defined as the FW–H acoustic source surface. These data are then used by the FW–H solver to compute the acoustic pressure at the specified receiver locations. Finally, the acoustic pressure signals are transformed into the frequency domain using the Fast Fourier Transform (FFT) to obtain the sound pressure level (SPL) spectra and the overall sound pressure level (OSPL).

2.5 Grid Independence and Validation Study

Due to the transient nature of the simulation, mesh independence and time-step sensitivity analyses are performed to assess the reliability of the numerical model. The effects of grid and time-step sizes are evaluated by comparing the energy extraction characteristics of three cases (D1–D3), as summarized in Table 2. The simulations are performed for the parameters: Ho/c = 1, θo = 70°, f* = 0.15, and Re = 8.58 × 104. The maximum relative error in the average power coefficient C P ¯ between cases D2 and D3 is approximately ~0.3%, indicating that further mesh refinement has a negligible effect on the numerical solution. Therefore, the D2 case, consisting of approximately 8 × 104 cells with a time step of 0.001T, is selected for all subsequent simulations to ensure an optimal balance between computational accuracy and cost.

Due to the lack of available experimental aeroacoustic data for the present operating conditions, validation is performed only for the aerodynamic and power extraction performance. Therefore, two validation studies are conducted to verify the accuracy and reliability of the present numerical methodology. The first validation is carried out for an arc-flapping motion case, as shown in Fig. 2a. The simulation is performed for the parameters: modified heaving amplitude Ho = 30°, radius-to-chord ratio R/c = 3, pitching amplitude θo = 50°, reduced frequency f* = 0.02, and Reynolds number Re = 5 × 105. The obtained results are compared with the existing data reported by Ma et al. [16]. As shown in Fig. 2b, the force coefficients (CV and CH) predicted in the present study closely match the results of Ma et al. [16]. The comparison demonstrates that the present numerical model is capable of accurately capturing the unsteady aerodynamic behavior of the oscillating foil.

A second validation study is performed by comparing the numerical results of the oscillating foil undergoing combined heaving and pitching motions with the experimental results of Simpson et al. [36] and the numerical results of Xu et al. [37]. The simulation is conducted for the following parameters: Ho/c = 1.23, θo = 85.82°, f* = 0.16, and Re = 1.38 × 104. As shown in Fig. 2c, the time history of lift coefficient obtained from the present simulation shows a good agreement with both the experimental data of Simpson et al. [36] and the numerical results of Xu et al. [37]. The maximum deviation between the present numerical results and the experimental data is approximately 10%. This demonstrates that the numerical model adopted in the present study is valid and reliable.

Figure 2: Validation of the present numerical model against existing numerical data (Ma et al. [16] and Xu et al. [37]) and experimental data (Simpson et al. [36]).

2.6 TOPSIS Multi-Criteria Decision Analysis

To identify the optimal operating condition by considering both energy harvesting performance and aerodynamic noise, the TOPSIS was employed. TOPSIS provides a systematic approach to rank the operating conditions by accounting for the trade-off between these two conflicting objectives. In the present study, C P ¯ was considered as the benefit criterion, whereas the maximum OSPL measured at receiver 5, where the highest acoustic levels were consistently observed, was considered as the cost criterion. Receiver 5 was selected because it consistently recorded the highest OSPL among all receiver locations, thereby representing the worst-case acoustic performance. Equal weighting (w = 0.5) was assigned to both criteria, giving equal importance to energy harvesting performance and aerodynamic noise. The TOPSIS method adopted in this study is outlined below. A decision problem is considered with m alternatives (A1, A2, ..., Am) and n evaluative criteria (C1, C2, ..., Cn) [38,39].

• 1.The decision matrix (X) is defined as:

The element ???????????? indicates the performance value of the i-th alternative with respect to the j-th criterion as follows [40]:

• 2.To eliminate the influence of different units and scales among the evaluation criteria, the decision matrix was normalized using the vector normalization method, as given by [41]:

The resulting normalized matrix is given by R = [rij]m×n.

• 3.The normalized matrix is multiplied by the criterion weights (wj) to obtain the equation below [38]:

• 4.Let JB and JC represent the sets of beneficial and non-beneficial criteria, respectively. The ideal (best) solution is defined as [41]:

The negative ideal (worst) solution is defined as:

• 5.The Euclidean distance of each alternative from the positive and negative ideal solutions, as given below [38,39,40,41,42]:

• 6.The relative closeness Ci* to the ideal solution is computed from the below form [42]:

Finally, the operating conditions are ranked in descending order according to C i * , and the operating condition with the highest value is considered the optimal solution.

3 Results and Discussions

3.1 Variation of Average Power Coefficient ( C P ¯ )

The effects of θo and f* on the average power coefficient C P ¯ of the OWPE are illustrated in Fig. 3. In the present study, four different pitching amplitudes (θo = 70°, 75°, 80° and 85°) are considered over a wide range of reduced frequencies (0.08 ≤ f* ≤ 0.24). It is evident from Fig. 3 that increasing the pitching amplitude gradually enhances the average power coefficient up to θo = 80°, beyond which the improvement becomes marginal and C P ¯ remains nearly constant. This behavior indicates that excessively large pitching amplitudes do not contribute significantly to further power enhancement. For θo = 70° and 75°, increasing the reduced frequency results in a continuous increase in C P ¯ up to f* = 0.15, after which the power coefficient decreases. In contrast, for larger pitching amplitudes of θo = 80° and 85°, the maximum C P ¯ is attained at a slightly higher reduced frequency of f* = 0.17, followed by a decline at higher frequencies. The reduction in power coefficient beyond the optimum frequency can be attributed to the occurrence of strong flow separation. Overall, the results suggest that the combined influence of pitching amplitude and reduced frequency plays a crucial role in determining the power extraction performance of the OWPE.

Figure 3: Variation of average power coefficient C P ¯  with reduced frequency f*.

3.2 Variation of Overall Sound Pressure Level (OSPL)

The effects of θo and f* on the overall sound pressure level (OSPL) predicted at the near-field receiver locations of the OWPE are illustrated in Fig. 4. It can be observed from Fig. 4a–d that increasing the pitching amplitude θo has only a marginal influence on the OSPL values at the selected receivers. In contrast, the reduced frequency f* has a more significant effect on the predicted acoustic levels at these near-field receiver locations. As f* increases, the OSPL decreases slightly up to f* = 0.12, beyond which it increases continuously over the remaining frequency range considered in this study. This increase in OSPL at higher reduced frequencies may be attributed to the intensified unsteady flow structures and stronger vortices generated during rapid oscillatory motion.

Figure 4: OSPL measured at 20 receivers for different f*. (a), (b), (c) and (d) denote θ_o = 70°, 75°, 80° and 85°, respectively.

The OSPL values recorded at the receivers are found to vary within the range 52.7 dB ≤ OSPL ≤ 110.2 dB. It should be noted that these values correspond to acoustic predictions obtained at receiver locations positioned within the near-field region, with receiver distances ranging from 1c to 10c from the OWPE. Therefore, the reported OSPL represent acoustic levels at the specified near-field receiver locations and should not be interpreted as far-field community noise levels. Furthermore, it is evident from Fig. 4 that increasing the streamwise distance of the receiver from the OWPE results in a continuous reduction in OSPL, as observed from the receiver pairs (1–4, 5–8, 9–12, 13–16, and 17–20). This reduction reflects the attenuation of the predicted acoustic pressure with increasing receiver distance. By considering both the average power coefficient C P ¯ and the corresponding OSPL trends shown in Fig. 3 and Fig. 4, high-performance operating conditions are generally observed around θo = 75°−80° and f* = 0.15–0.17. This operating condition provides a favorable balance between efficient power extraction and reduced acoustic levels at the selected near-field receivers.

Although the present acoustic predictions are obtained using the FW–H acoustic analogy, the receiver locations considered in this study are located within the near-field region (1c–10c) surrounding the OWPE. Therefore, the reported OSPL values should be interpreted as predicted acoustic levels at the specified receiver locations rather than direct estimates of far-field community noise. For engineering applications, a first-order estimate of the far-field sound pressure level may be obtained by assuming spherical wave propagation, according to Lp,2=Lp,1−20log10r2r1(35) where Lp,1 is the SPL at the reference receiver distance r1, and Lp,2 is the estimated SPL at a farther receiver distance r2. This approximation accounts only for the effect of distance on sound propagation. A more comprehensive far-field analysis would also consider source directivity and other propagation effects.

3.3 Time Histories of Power Coefficient and Vortex Structures

The time histories of CY, CP, CPY and CPθ along with the corresponding vortex structures are compared in Fig. 5 for two different operating conditions. The first case corresponds to θo = 75° and f* = 0.15, while the second case corresponds to θo = 80° and f* = 0.17. Each time period comprises of two strokes, downstroke (t/T = 0−0.5) and upstroke (t/T = 0.5−1). The instantaneous power coefficient (CP) exhibits two local maxima and two local minima during each time period. The improvement in CP primarily noticed in the middle of downstroke and upstroke. This increase in CP is primarily attributed to the development and evolution of the leading-edge vortex (LEV), which plays a crucial role in the dynamic stall mechanism. To comprehensively analyze the flow evolution during one oscillation cycle, the instantaneous z-vorticity and streamwise velocity contours for the first case (θo = 75° and f* = 0.15) are presented in Fig. 6 and Fig. 7, respectively.

Figure 5: Time histories of force and power coefficients (a) θ_o = 75° and f* = 0.15 (b) θ_o = 80° and f* = 0.17.

Figure 6: Vorticity with pressure isolines. (a) t/T = 0 (b) t/T = 0.1 (c) t/T = 0.2 (d) t/T = 0.3 (e) t/T = 0.4 (f) t/T = 0.5 (g) t/T = 0.6 (h) t/T = 0.7 (i) t/T = 0.8 (j) t/T = 0.9 (k) t/T = 1.

At the beginning of the oscillation cycle (t/T = 0), a separated vortex can be observed near the trailing edge (TE) on the upper surface of the wing (Fig. 6). As the wing starts its downstroke, a new leading-edge vortex (LEV) begins to develop on the lower surface of the wing, which can be clearly seen at time instant t/T = 0.2. The LEV continues to strengthen in the time interval t/T = 0−0.3 as the wing moves down. During this interval, the LEV remains attached to the wing surface, creating a strong low-pressure region that enhances the lift and consequently increases the power extraction, as observed in Fig. 5. At the end of the downstroke (t/T = 0.4−0.5), the LEV begins to detach from the wing surface. As the wing reverses from the downstroke to the upstroke, a new LEV with opposite vorticity begins to develop on the upper surface. Similar flow features are observed during the upstroke, where the LEV grows, convects downstream along the wing surface, and eventually sheds into the wake.

The streamwise velocity contours presented in Fig. 7 further illustrate the flow evolution during the oscillation cycle and provide more insight into the development and shedding of the LEV. As the LEV develops, a region of accelerated flow appears near the leading edge. During the downstroke, the LEV continues to grow and strengthen, resulting in a stronger velocity gradient between the upper and lower surfaces of the wing. After the LEV detaches, the wake becomes wider and a larger low-velocity region develops behind the wing. This wake structure results from the interaction between the detached LEV and the trailing-edge vortex shed during the previous cycle.

Figure 7: Velocity contour. (a) t/T = 0 (b) t/T = 0.1 (c) t/T = 0.2 (d) t/T = 0.3 (e) t/T = 0.4 (f) t/T = 0.5.

As discussed earlier, a comparison of the two operating conditions (θo = 75° and f* = 0.15 and θo = 80° and f* = 0.17), as shown in Fig. 5, indicates that the instantaneous power coefficient (CP) is higher for the second operating condition. Fig. 6 and Fig. 7 demonstrate that the formation, growth, and shedding of the LEV strongly influence the instantaneous aerodynamic forces and the power generation characteristics of the OWPE. A comparison of the flow fields for the two operating conditions in Fig. 8 shows that the second operating condition generates stronger vortex structures around the wing and in the wake. The stronger LEV enhances the suction region over the wing surface, resulting in increased aerodynamic lift and consequently a higher instantaneous power coefficient (CP) than the first operating condition. However, the stronger vortices also generate larger unsteady pressure fluctuations around the wing and in the wake. These pressure fluctuations contribute to the predicted acoustic response at the near-field receiver locations, leading to higher OSPL values at the selected receivers.

Figure 8: Instantaneous vortex structures (z-vorticity) at the selected time instants indicated by the colored markers in the time-history plots shown in Fig. 5 (a) θ_o = 75° and f* = 0.15 (b) θ_o = 80° and f* = 0.17.

In addition to the differences in vortex strength, an increase in flapping frequency also affects the aeroacoustic characteristics predicted at the receiver locations. A higher flapping frequency alters the generation and shedding of vortices. Consequently, the pressure fluctuations around the wing increase, resulting in higher predicted acoustic levels at the near-field receivers. Although the development of the LEV enhances the lift force generation and contributes to the increase in instantaneous power coefficient, the higher OSPL values predicted at the selected receiver locations are partly attributed to the increased flapping frequency and the associated unsteady pressure fluctuations.

3.4 Variation of Sound Pressure Level (SPL) with Acoustic Frequency (F_a)

The spectral analysis of the predicted acoustic pressure at the receiver locations, represented in terms of sound pressure level (SPL) over a wide range of acoustic frequencies (fa), are illustrated in Fig. 9. In the present study, the maximum acoustic frequency considered is 2400 Hz. The results demonstrate that the reduced frequency f* significantly influences the SPL distribution at the receiver locations. As the reduced frequency f* increases, the SPL spectra maintain relatively higher SPL values over a broader range of acoustic frequencies. For example, comparing f* = 0.08 and f* = 0.15, the case with f* = 0.24 exhibits consistently higher SPL values over a wider range of acoustic frequencies, thereby increasing the corresponding overall sound pressure level (OSPL). In several receiver groups, especially for f* = 0.15 and f* = 0.24, a broadband hump appears in the medium-frequency range, suggesting that the acoustic response is not purely tonal but contains distributed broadband components. Furthermore, for a given reduced frequency, the SPL generally decreases with increasing receiver distance from the OWPE, consistent with the attenuation of the acoustic pressure as the sound propagates away from the source.

Irrespective of f* and receiver location, the SPL is generally higher at lower acoustic frequencies and decreases progressively with increasing acoustic frequency. This indicates that the predicted acoustic response at the near-field receiver locations of the OWPE is mainly dominated by low-frequency components. Although broadband fluctuations are present at higher acoustic frequencies, the SPL remains bounded and does not exceed the dominant low-frequency levels. Furthermore, the SPL measured at all defined receivers remains below 55 dB for fa >1000 Hz, indicating that the high-frequency contribution is relatively weak compared with the low-frequency components.

Figure 9: SPL (dB) measured at 20 receiver locations for different f*. (a) R1−R4 (b) R5−R8 (c) R9−R12 (d) R13−R16 (e) R17−R20.

3.5 Identification of the Optimal Operating Condition Using the TOPSIS Method

Table 3 summarizes the TOPSIS ranking for all 52 operating conditions. The results indicate that the operating condition with θo = 80° and f* = 0.16 the highest TOPSIS score of 0.9553, followed by θo = 80° and f* = 0.17 and θo = 80° and f* = 0.15, with TOPSIS scores of 0.9443 and 0.9433, respectively. These operating conditions provide the best compromise between maximizing the average power coefficient and minimizing the aerodynamic noise. The TOPSIS ranking further shows that the highest-ranked operating conditions are predominantly concentrated around θo = 80° and f* = 0.15–0.18, indicating that this operating range offers the most favorable balance between the two objectives.

4 Conclusion

The present numerical study systematically investigates the effects of the pitching amplitude θo and reduced frequency f* on the power extraction performance and aeroacoustic characteristics of an OWPE. Numerical simulations are performed for a NACA 0015 airfoil undergoing prescribed heaving and pitching motions, and the unsteady flow field is resolved using a dynamic mesh with smoothing and local remeshing. The aerodynamic noise is predicted using the Ffowcs Williams–Hawkings (FW–H) acoustic analogy model at 20 receiver locations positioned in the near-field region surrounding the OWPE. The results show that the reduced frequency f* has a strong influence on the near-field overall sound pressure level (OSPL), while the effect of the pitching amplitude θo is relatively less significant. Increasing the reduced frequency broadens the acoustic spectrum and increases the predicted near-field sound pressure, while the acoustic response remains dominated by low-frequency components. The predicted SPL also decreases with increasing receiver distance due to acoustic attenuation. The reported results are limited to the near field while the far-field acoustic levels may be estimated using a spherical wave propagation model. The flow-field analysis further reveals that the formation, growth, and shedding of the leading-edge vortex (LEV) strongly influence both the power extraction and the near-field aeroacoustic characteristics of the OWPE. Stronger LEVs enhance power extraction but also increase unsteady pressure fluctuations, resulting in higher near-field acoustic levels. Finally, the TOPSIS-based multi-criteria analysis identified θo = 80° and f* = 0.16 as the optimal configuration, providing the best compromise between power extraction and near-field aerodynamic noise. The findings of this study provide a better understanding of the relationship between unsteady vortex dynamics, aerodynamic performance, and aeroacoustic characteristics of OWPEs. These results provide useful guidance for designing quieter and more efficient OWPEs for urban wind energy applications.