Constructive interference in a network of elastically-bounded flapping plates

TL;DR

This study investigates how arrays of elastically-bounded flapping plates can cooperatively enhance energy harvesting efficiency through constructive interference, combining numerical simulations and wind-tunnel experiments.

Contribution

It introduces a comprehensive analysis of multi-device aeroelastic energy harvesters, revealing beneficial cooperative effects and optimal configurations for power amplification.

Findings

Constructive interference can double power output in side-by-side arrangements.

Elastic tuning improves downstream device performance.

Numerical results are validated with wind-tunnel experiments.

Abstract

Aeroelastic phenomena are gaining significant attention from the perspective of energy harvesting (EH) with promising applications in supplying low-power remote sensors. Besides the development of individual EH devices, further issues are posed when considering multiple objects for realizing arrays of devices and magnifying the extracted power. Due to nonlinear mutual interactions, the resulting dynamics is generally different from that of single devices and the setup optimisation turns out to be nontrivial. In this work, we investigate the problem focusing on a flutter-based EH system consisting of a rigid plate anchored by elastic elements and invested by a uniform laminar flow, undergoing regular limit-cycle oscillations and flapping motions of finite amplitude. We consider a simplified, yet general, physical model and employ three-dimensional direct numerical simulations based on a finite-difference Navier-Stokes solver combined with a moving-least-squares immersed boundary method. Focusing on main kinematic and performance-related quantities, we first report on the dynamics of the single device and then on multiple devices. A parametric exploration is performed by varying the mutual distance between the devices. For the in-line arrangement, a recovery in performance for downstream devices is achieved by tuning their elasticity. Moreover, cooperative effects in the side-by-side arrangement are found to be substantially beneficial in terms of resulting power, with increases (i.e. constructive interference) up to 100% with respect to the single-device configuration. In order to confirm this numerical evidence, complementary results from wind-tunnel experiments are presented. Finally, we describe the system behaviour when increasing further the number of devices, outlining the ultimate goal of developing a high-performance EH network of numerous aeroelastic energy harvesters.

Tables

Table 1. Table 1: Performance of side-by-side arrays with different number of devices N d subscript 𝑁 d N_{\mathrm{d}} ; C ¯ P tot superscript subscript ¯ 𝐶 𝑃 tot \bar{C}_{P}^{\mathrm{tot}} is the total power coefficient of the array (Eq. ( 10 )) and I ​ F 𝐼 𝐹 IF is the interference factor as defined by Eq. ( 11 ).

$N_d$	${\overline{C}}_P^{\mathrm{tot}}$	$IF$
3	1.10	0.76
5	2.01	0.92
11	4.56	0.98

Table 2. Table 2: Flapping observables for the experimental array of three devices in side-by-side arrangement. The suffix denotes values for the single (S) or the array (A) configuration.

Device	$\mathrm{Amp}{\left(y_{\mathrm{𝑃𝑃}}\right)}^S$	$\mathrm{Amp}{\left(y_{\mathrm{𝑃𝑃}}\right)}^A$	${\mathrm{𝑆𝑡}}^S$	${\mathrm{𝑆𝑡}}^A$	${\overline{C}}_P^S$	${\overline{C}}_P^A$
W2’	0.69	0.97	0.16	0.17	0.13	0.30
W1	0.73	1.07	0.16	0.17	0.16	0.43
W2	0.83	0.9	0.17	0.17	0.23	0.29

Table 3. Table 3: Grid settings used for convergence study. N x subscript 𝑁 𝑥 N_{x} , N y subscript 𝑁 𝑦 N_{y} and N z subscript 𝑁 𝑧 N_{z} denote the number of nodes in the x 𝑥 x , y 𝑦 y and z 𝑧 z direction, while H 𝐻 H and h ℎ h indicate the minimum and maximum resolution, respectively.

Grid	$N_x$	$N_y$	$N_z$	$H$	$h$
A	251	84	127	$0.4$	$0.04$
B	271	97	147	$0.35$	$0.035$
C	334	113	170	$0.3$	$0.03$
D	401	136	205	$0.25$	$0.025$
E	501	170	257	$0.2$	$0.02$
F	541	189	286	$0.18$	$0.018$
G	601	213	322	$0.16$	$0.016$
H	641	235	353	$0.146$	$0.0146$
I	751	258	388	$0.133$	$0.0133$
J	865	298	448	$0.115$	$0.015$
K	1001	344	516	$0.1$	$0.01$

Equations

$$m{\ddot{\mathbf{r}}}_G={\mathbf{F}}_{\mathrm{aero}}+{\mathbf{F}}_{\mathrm{el}},$$

$$I_x\ddot{\theta }=M_{\mathrm{aero}}+M_{\mathrm{el}},$$

$${\mathbf{F}}_{\mathrm{aero}}={\int }_S(\mathbf{\tau }-p\mathbf{n})\mathrm{d}S,$$

$$M_{\mathrm{aero}}={\int }_S\mathbf{r}\times (\mathbf{\tau }-p\mathbf{n})\mathrm{d}S,$$

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$$K_{\mathrm{cr}}\approx \frac{3}{4}\pi \text{AR}.$$

$${\bar{C}}_P=\frac{{\bar{P}}_y}{\frac{1}{2}{\rho }_{\mathrm{f}}sc{U}^3},$$

$$\eta =\frac{{\bar{P}}_y}{\frac{1}{2}{\rho }_{\mathrm{f}}sd{U}^3},$$

$${\bar{C}}_P^{\mathrm{tot}}=\sum _{i=1}^{N_{\mathrm{d}}}{\bar{C}}_P=N_{\mathrm{d}}{\bar{C}}_P^{\mathrm{S}}+I,$$

$$IF=\frac{I}{N_{\mathrm{d}}{\bar{C}}_P^{\mathrm{S}}}.$$

$$f(h)\approx f^{\ast }+C{h}^p$$

$$ϵ(h)=\frac{f(h)-f(h_{\min })}{f(h_{\min })},$$

$$ϵ(h)=C(h^p-h_{\min }^p),$$

Full text

††journal: Journal of Fluid and Structures

1 Introduction

Flow-induced vibrations are receiving growing interest from the perspective of energy harvesting (EH), in view of applications related to the supply of low-power sensor networks within the framework of the Internet of Things [1, 2]. A variety of novel devices and corresponding models has been proposed in literature and they can be classified depending on the aeroelastic instability occurring, such as, e.g., flutter, galloping and vortex-induced vibrations. In all cases, the understanding of these nonlinear phenomena is crucial in order to design efficient and reliable EH devices and to fine-tune their performance.

For slender aerodynamic bodies such as wings or plates, the typical mechanism to be exploited is aeroelastic flutter, consisting of a self-excitation between flow and structural response (see, e.g., Ref. [2] for an introduction to flutter-based EH). Among the various flutter-based EH systems, one can include those based on flapping flags [3, 4, 5] as well as passively flapping airfoils [6, 7]. For the latter, the motion of an essentially rigid, streamlined body is limited to a certain number of degrees of freedom (DoFs), typically two: translation along the transverse direction (plunge) and rotation around a spanwise axis (pitch).

Focusing on fully-passive systems (i.e., those whose dynamics is entirely governed by the fluid-structure interaction, without any prescribed kinematics), early research by Peng and Zhu [8], conducted by two-dimensional numerical simulations, highlighted the variety of possible flapping states and provided a first estimate of the resulting power and efficiency. These findings were confirmed and enriched by further computational studies in closely comparable conditions [9, 10] as well as rather different systems, e.g. with different cross sectional shapes and/or operating at higher-Reynolds regimes [11, 12, 13, 14].

Young et al. [11] introduced a configuration where the pitch and plunge motions are constrained by a mechanical linkage along with pitch control in order to increase the performance. Veilleux and Dumas [12] performed an optimization study for a fully-passive device by two-dimensional CFD simulations which became the basis for the experimental prototype later presented by Boudreau et al. [15]. These studies agree in reporting that the performance is improved for an adequate synchronization between the two DoFs and in case of nonsinusoidal pitching motion. Wind-tunnel investigations were reported by Pigolotti et al. [16, 17, 18] considering a flat plate in a classical pitch-and-plunge arrangement and exploring systematically the effect of several physical parameters on the flutter onset and the nonlinear oscillations in the postcritical regime. A nearly identical system was considered in the work by Wang et al. [14] where two-dimensional computations were performed in order to characterize the dependence of limit-cycle oscillations with respect to the governing parameters, initial conditions, spring nonlinearity and extraction (modelled by viscous damping).

Besides the development of the individual energy harvester and the consequent understanding of the associated nonlinear system, a further degree of complexity arises when looking at a network of multiple energy harvesters, that allow to scale the total extracted power. In fact, the behaviour of each device of the network is coupled to all the others due to the nonlinear character of the dynamical system.

Among the few contributions given in this latter perspective, Bryant et al. [19] tested devices made by a rigid airfoil hinged on a cantilevered flexible beam and arranged in in-line or staggered configuration, reporting performance improvements for downstream devices due to the beneficial effect of wake forcing. Moreover, it was later proposed that such effect can be controlled by tuning the pitching stiffness of downstream devices [20]. McCarthy et al. [21, 22] focused on the in-line arrangement considering a different flutter-based system where a triangular leaf is joined to a piezoelectric stalk by a revolute hinge. They reported increases in power of about $40\%$ for the downstream device caused by the interaction with the horseshoe cone vortex released by the upstream device.

Even if they have a different behaviour, flapping flags were also investigated in the same framework. Several works focused, in particular, on multiple filaments or flexible plates placed in the side-by-side configuration, showing increased oscillations over an intermediate range of distances between the flags [23, 24, 25].

In this work, we focus on an aeroelastic system based on elastically-anchored plates experiencing a fluttering instability when invested by laminar flow, giving rise to finite-amplitude limit cycle oscillations (LCOs) and flapping motions [26, 27, 28, 29, 30]. Previous related work addressed the identification of the critical threshold for sustained flapping [27, 28] and a first characterization of the postcritical behaviour for a representative configuration of the real device [29], as well as the development of an experimental prototype equipped for energy extraction by electromagnetic coupling [30]. Overall, we aim at developing EH devices of centimetric size able to extract electrical power $\mathcal{O}(\mathrm{mW})$ in low wind conditions (i.e., $U<5\,\mathrm{m/s}$). Despite the analogies with flutter-based systems, differences exist with respect to the classical pitch-and-plunge model, the main being the absence of rotational elastic constraints yielding a different dynamics, e.g. in the critical condition for flapping onset [28].

The goal of the present work is to provide a fundamental study on the behaviour of multiple elastically-bounded flapping plates immersed in an incompressible laminar flow, characterizing the resulting dynamics and giving useful insights for the development of EH networks made of arrays of such devices. To this aim, we will employ three-dimensional numerical simulations combining a finite-difference Navier-Stokes solver and an immersed boundary technique. First we will characterize the dynamics of the single device, and then investigate three basic arrangements of multiple devices, i.e.: (i) in-line, (ii) staggered and (iii) side-by-side. A parametric study will be conducted by focusing on the dependence of main quantities of interest with respect to the mutual distance between devices. In order to corroborate the numerical evidence, we will provide experimental results from wind-tunnel measurements pertaining to real EH applications.

Following this introduction, the rest of the paper is structured as follows. In Sec. 2 we introduce the representative aeroelastic system and the main governing parameters; Sec. 3 concerns the numerical method, while results are provided in Sec. 4. Conclusions are drawn in Sec. 5. In a final Appendix we provide evidence for the validation of the numerical method and for all the run parameters used for the simulations.

2 Physical model

We consider a rigid plate of homogeneous density $\rho_{\mathrm{s}}$ and geometry characterized by chord $c$, span $s$ and thickness $\delta$. The plate is restrained by two linear springs with stiffness $k$ and zero restlength. As shown in Fig. 1a, for each spring one end is connected to the plate at the pivot point (PP), moving with the body, while the other one is retained fixed at the anchor point (AP). The body is invested by a uniform flow with inflow velocity $U$, the fluid has kinematic viscosity $\nu$ and density $\rho_{\mathrm{f}}$. As additional assumptions, gravity effects are neglected and energy extraction is not considered since we are concerned with the purely aeroelastic behaviour of the system.

Looking at the introduced quantities, the following four nondimensional parameters can be derived: the chord-based Reynolds number $\mathit{Re}=Uc/\nu$, the nondimensional stiffness $K=2\,k/(\rho_{\mathrm{f}}U^{2}c)$ (accounting for both springs), the two-dimensional mass parameter $\rho_{\mathrm{w}}=\rho_{\mathrm{s}}\delta/(\rho_{\mathrm{f}}c)$ (where the subscript ‘w’ stands for wing, this quantity depending both on the density ratio and the wing cross section) and the planform aspect ratio $\AR=s/c$. In the rest of the work we will always refer implicitly to nondimensional quantities, including spatial and temporal coordinates, i.e.: $(x,y,z)\mapsto(x/c,y/c,z/c)$, $t\mapsto t/(c/U)$.

Three degrees of freedom (DOFs) are allowed in the plate motion: (i) translation in the streamwise direction (surge), (ii) translation in the transverse direction (plunge), (iii) rotation around the spanwise axis (pitch). Under these assumptions, the rigid body equations governing the motion of the center of mass (whose position is denoted by $\mathbf{r}_{G}$) and rotation $\theta$ (written with respect to $G$) read:

[TABLE]

[TABLE]

where $m=\rho_{\mathrm{w}}\,\AR$ is the mass, $I_{x}$ is the moment of inertia with respect to an axis passing through the hinge and directed in the spanwise direction, $\mathbf{F}_{\mathrm{aero}}$ is the aerodynamic force, $\mathbf{F}_{\mathrm{el}}=-K\,\mathbf{r}_{PP}$ is the elastic force exerted by the springs, ${M}_{\mathrm{aero}}$ is the aerodynamic moment and ${M}_{\mathrm{el}}=(\mathbf{r}_{PP}-\mathbf{r}_{G})\times\mathbf{F}_{\mathrm{el}}$ is the elastic one. The aerodynamic force and moment are obtained by integrating the fluid stress tensor over the plate surface $S$:

[TABLE]

[TABLE]

where $\boldsymbol{\tau}$ is the viscous shear stress, $p$ is the pressure, $\mathbf{n}$ is the unit vector normal to the plate and $\mathbf{r}$ is the distance with respect to the center of mass $G$.

Although differences exist with respect to the real device (e.g., in mass distribution), this kind of description, already considered in [27, 28] and extended here to the three-dimensional case, is able to reproduce the essential physics that has been observed experimentally and thus matches our goal of understanding, at least qualitatively, the peculiar effects arising when multiple devices are considered.

3 Numerical method

The flow obeys to the incompressible Navier-Stokes equations which, in nondimensional form, read:

[TABLE]

where $\mathbf{u}=\mathbf{u}(\mathbf{x},t)$ is the fluid velocity, $p=p(\mathbf{x},t)$ the pressure and $\mathbf{f}=\mathbf{f}(\mathbf{x},t)$ a volumetric forcing. We consider a three-dimensional domain of size $L_{x}\times L_{y}\times L_{z}$ (Fig. 1b), with the following boundary conditions: the fluid velocity is uniform at the inlet, convective boundary conditions are used at the outlet [31], the top and bottom faces are treated as slip (i.e., non-penetrating) walls, while periodicity is assumed at side faces.

Since we aim at dealing with multiple moving objects, we resort to the immersed boundary (IB) technique [32] and employ the moving-least-square (MLS) method recently proposed by de Tullio and Pascazio [33] that here we briefly summarize (for more details, see Refs. [34, 33, 35]).

Eqs. (5) and (6) are solved numerically on a Cartesian grid, with the forcing term $\mathbf{f}$ mimicking the presence of solid bodies. Centered finite differencing is adopted for space discretization using a staggered grid, the overall scheme yielding second-order accuracy. The main iteration loop can be summarized as follows: First, the nonlinear terms are computed explicitly using a second-order Adams-Bashfort scheme. Next, the diffusive terms are discretized implicitly with a Crank-Nicolson scheme and an intermediate, nonsolenoidal velocity field is computed by an approximate factorization technique for the resulting algebraic system [34]. The correction by the IB forcing is then added to the velocity. Finally, the Poisson equation for mass conservation is solved and the divergence-free velocity is obtained, together with the pressure field.

Concerning the IB treatment, a moving-least-square (MLS) interpolation is used to reconstruct the solution at the immersed surface [36]. The plate is discretized by a planar surface mesh with $N_{\mathrm{t}}$ triangular elements, their centroids being the Lagrangian markers at which the forcing is computed by imposing the no-slip condition. For each Lagrangian marker, we consider a support domain enclosing $N_{\mathrm{e}}=27$ adjacent Eulerian nodes which are used for the interpolation and spreading operations. Compared to classical IB approaches, the adoption of the described method allows to use a larger simulation timestep and to obtain smoother flow solutions and aerodynamic loads.

The numerical procedure has been tested for the present application by a convergence analysis with respect to the spatial and temporal resolution. Results of this latter and discussion on the choice of grid parameters are collected in A. As a baseline indication, the used domain box ranges from $(-5,-5,-5)$ to $(5,5,10)$ and is discretized using a minimum grid spacing $h=0.02$, with the triangulated mesh discretizing the solid plate having comparable resolution, so that for the case $\AR=2$ the total number of triangular elements is $N_{\mathrm{t}}\approx 1.3\times 10^{4}$.

4 Results

4.1 Single device

As a starting point, we characterize the dynamics of the single, isolated device. Two points are investigated to this aim: (i) the identification of the critical condition for sustained flapping and (ii) the dependence of the limit-cycle features on the governing parameters. The first point can be addressed by recalling the predictive arguments of Ref. [28], based on the coupling between the natural frequency and the frequency of the pitching response for small angles in the case of hinged plate. Recasting these arguments for our model, the following expression is found for the critical $K$:

[TABLE]

In order to verify this prediction, we have performed a series of simulations at $\rho_{\mathrm{w}}=2$ and $\mathit{Re}=100$ while varying $\AR$ and $K$, whose results are collected in Fig. 2. The plot shows that indeed sustained flapping occurs for $K<K_{\mathrm{cr}}$ (unstable cases); on the other hand, for $K>K_{\mathrm{cr}}$ the wing asymptotically aligns with the flow (stable cases). As predicted by the theory, a linear dependence of $K_{\mathrm{cr}}$ on $\AR$ is found. Although here the simulations have been performed at $\mathit{Re}=100$, the same evidence was found for different values (tested up to $\mathit{Re}=1000$), in agreement with Eq. (7) where the Reynolds number does not explicitly appear. Furthermore, similar behaviour has been verified when varying the density parameter $\rho_{\mathrm{w}}$.

Let us now move on to the analysis of the post-critical, nonlinear regime, fixing the aspect ratio to $\AR=2$ and exploring the dependence of the flapping dynamics for different values of the nondimensional stiffness ($K=\{1,2,3,4\}$) and Reynolds number ($\mathit{Re}=\{100,200,500,1000\}$).

The resulting system kinematics and performance are summarized in Fig. 3. The first quantity of interest (also from the EH viewpoint) is the peak-to-peak amplitude of the PP displacement along the transverse distance. As shown in Fig. 3a, this is found to increase for decreasing $K$ (i.e., softening the spring). On the other hand, the Strouhal number (defined as $\mathit{St}=fc/U$, where $f$ is the flapping frequency) increases almost linearly with $K$ (Fig. 3b).

Concerning the performance parameters, we refer to the average plunge power coefficient $\bar{C}_{P}$ and the Betz efficiency $\eta$, defined as:

[TABLE]

[TABLE]

where $\bar{P}_{y}=\frac{1}{T}\int_{T}F_{\mathrm{aero}}^{y}\,\dot{y}_{PP}\,\mathrm{d}t$ is the average power associated with the plunge motion, $F_{\mathrm{aero}}^{y}$ is the vertical component of the aerodynamic force (i.e., the lift force), and $d$ is the maximum transverse distance swept by the wing during its motion. Looking at Fig. 3c, the trend of the power coefficient appears to be qualitatively similar to those found for the PP amplitude. Conversely, the Betz efficiency exhibits a trend which is not monotonic (Fig. 3d).

Overall, variations of the Reynolds number do not appear to modify the described trend. For lower $\mathit{Re}$, however, the oscillation amplitude (as well as power and efficiency) weakens, consistently with a more dominant effect of viscosity, except for the stiffest case $K=4$ where for $\mathit{Re}=1000$ we observe a sharp decrease. In this case, it can be observed that the wake remains attached to the wing, while for lower $\mathit{Re}$ a sequence of vortices is released from the leading edge, producing higher oscillations.

In light of this evidence, in the following we will focus uniquely on the case at $\mathit{Re}=100$. This value is selected to deal with the smoothest flow solution among the considered cases, in order to get a clearer understanding of the basic mechanisms occuring in the interaction between multiple devices.

4.2 In-line arrangement

We begin our study on the interaction of multiple devices by considering an in-line configuration, i.e. a second wing is placed downstream at a distance $\mathbf{r}=(0,0,r_{z})$ along the streamwise direction, as sketched in Fig. 4a. We retain the initial perturbation only for the upstream wing to evidence how the downstream wing dynamics is affected by the impacting wake. To accommodate a second downstream device, the domain is enlarged in the streamwise direction up to $z=20$.

To start our analysis, we fix $K=3$ for both devices, since this was found to be the most efficient condition (i.e., maximising $\eta$) in the case of the single device, and perform simulations for different values of the relative distance $r_{z}$ (see Fig. 4a). An insight of the resulting dynamics for the case $r_{z}=2$ is given by the top panels of Fig. 5. Looking at the time trace of the plunging motion (Fig. 5a), after a short transient (about 2 cycles) a phase shift is established between the oscillations of the two devices. In the new flapping state, the motion of the downstream wing gets synchronized to the wake released by the upstream one, the pitching motion being driven by low-pressure vortical regions, as it can be observed from Fig. 6, showing instantaneous views of the vorticity field within one flapping cycle (see also Movie 1 provided with the Supplementary Material). Furthermore, from the shape of limit-cycles reported in Fig. 5b, we note that the state-space trajectory of the second device is contained within that of the upstream device.

Similar findings are obtained when varying the distance $r_{z}$, as shown in Fig. 7: the downstream device always oscillates with amplitude smaller than the upstream one (Fig. 7a), with a monotonic trend that seems to recover the single-device behaviour for large $r_{z}$, as expected. The flapping frequencies of the two devices, shown in Fig. 7b, are essentially locked and slightly smaller than that of the single case. $\mathit{St}$ decreases while decreasing the separation distance, reflecting an alteration of the upstream wing dynamics, as well. Finally, the power coefficient shows a similar trend (Fig. 7c), although the difference with respect to the single configuration is more pronounced: in fact, this quantity involves the product between the lift force and the PP velocity, both being weakened.

4.2.1 Structural tuning of downstream device

In order to improve the performance of the downstream device, we investigate the configuration in which the value of its stiffness $K_{(2)}$ is varied, while for the upstream device it is mantained fixed to $K_{(1)}=3$. Two streamwise distances are considered, i.e. $r_{z}=\{2,4\}$. Similarly to the previous case, the flapping observables are reported in Fig. 8. The plunge motion increases its amplitude for decreasing $K_{(2)}$, attaining values close to the device in the single configuration (Fig. 8a). Conversely, the Strouhal number (Fig. 8b) decreases while softening the spring for $K_{(2)}<K_{(1)}=3$. In this range, the resulting dynamics is a nonlinear combination of the self-excitation, taking place in uniform flow, and the wake-forcing mechanism previously discussed.

Bottom panels of Fig. 5 show the resulting flapping in time for a representative case, where one can notice the larger amplitude of both pitch and plunge compared to the corresponding case with $K_{2}=K_{1}$. Furthermore, for this particular case ($r_{z}=2$, $K_{2}=1$), one can see that the system shows asymmetric and multiperiodic oscillation of the pivot-point (Fig. 5c,d). Indeed, it was observed that when $K_{(2)}<K_{(1)}$ the spectral content is richer and the dominant frequency is altered compared to the single-wing case at $K=K_{(2)}$.

For $K_{(2)}>K_{(1)}=3$, the flapping frequency remains locked to the upstream wing. The downstream-wing oscillation decreases while increasing $K_{(2)}$ but flapping now occurs also for $K_{(2)}>K_{\mathrm{cr}}\approx 4.7$, unlike what happens in the case of an isolated device. In this case, the only mechanism causing such motion is wake forcing.

A recovery in the value of the power coefficient is obtained by lowering $K_{2}$, as shown in Fig. 8c, which is beneficial from the EH perspective, although in none of the considered cases the same amount as for the single device was obtained. Our numerical results can be compared with the experimental evidence recently presented in Ref. [20], where the concept of tuning the pitching stiffness of the downstream wing was proposed for pitch-and-plunge EH systems: despite the different structural features (e.g., the presence of restoring moments both in plunging and pitching) and the significant difference in Reynolds numbers considered, a qualitative analogy can be drawn.

4.3 Staggered arrangement

In this case (Fig. 4b), a pair of devices (W2 and W2’) is placed aft of the first one (W1) at a distance $\mathbf{r}=(0,\pm r_{y},r_{z})$. We choose the closest streamwise distance considered for the in-line case, i.e. $r_{z}=2$, and investigate the behaviour for different values of the transverse distance $r_{y}=\{0.5,1,1.5,2\}$.

Results of the analysis are presented in Fig. 9.

Concerning the plunging amplitude (Fig. 9a), an attenuation is found with respect to the single-device configuration for all devices, which is however limited for W1 to about $5\%$, while for W2 and W2’ we have a relative peak at $r_{y}=1$ (approaching the value in single configuration) and a sharp decrease for $r_{y}=0.5$ (about $20\%$). In the latter, the downstream wings are found to lie entirely within the wake released by W1, the oscillation being reduced in a way similar to what reported for the in-line arrangement, featuring the same synchronization mechanism. For $r_{y}=1$, the wake is impacting on downstream devices only during a portion of the flapping cycle, yielding a different dynamics: the flapping motions are collectively in-phase, as shown by Fig. 10, and the occurring interference is presumably responsible for the relative peak in amplitude. Increasing the distance to $r_{y}=1.5$ and 2, the same in-phase motion is still observed but W2 and W2’ are now located outside the wake region and consequently the weakly constructive interference is not present anymore.

Note that for W2 and W2’ the flapping motion is asymmetric with respect to the streamwise direction but substantially specular with respect to each other. The Strouhal number shows almost negligible variations (less than $5\%$) compared to the single-device configuration, decreasing for smaller $r_{y}$ (Fig. 9b). Values of the power coefficient are always found to be lower than that of the isolated device. We note however that, for $r_{y}>1$, the power coefficient is found to be higher for downstream devices (Fig. 9c).

4.4 Side-by-side arrangement

We now focus on the side-by-side arrangement, where we place a pair of devices (W2 and W2’) at a transverse distance $\mathbf{r}=(0,\pm r_{y},0)$ with respect to the central one (W1), as sketched in Fig. 4c. As for the other two arrangements, only the first wing is initially perturbed.

Flapping observables from numerical simulations considering three different values of the mutual distance, i.e. $r_{y}=\{1,1.5,2\}$, are reported in Fig. 11. Looking at the plunging amplitude (Fig. 11a), for sufficiently small separations (i.e. $r_{y}\leq 1.5$) the resulting values are larger than that obtained for the single device, although the relative increment is only up to about $8\%$ for W1 and $20\%$ for W2 and W2’. The flapping frequency is nearly identical for the three devices and slightly increases while decreasing $r_{y}$, up to about $5\%$ when $r_{y}=1$ (Fig. 11b). Despite the relatively small variation of these two quantities, a significant increase of the power coefficient occurs for all flapping plates (Fig. 11c): with respect to the single device, for the central wing the increase is almost $50\%$, while for the side wings this reaches nearly $90\%$.

Similarly to what it was observed for the staggered arrangement (Sec. 4.3), the central wing undergoes symmetrical motion while the side wings exhibit slightly asymmetrical flapping, as it can be observed from the time traces of the transverse PP displacement reported in Fig. 12 (see also Movie 2 in the Supplementary Material). Moreover, the steady-state oscillations of the side wings are approximately in-phase with respect to each other and in counter-phase with respect to the central wing. From the qualitative viewpoint, we can observe how the resulting scenario resembles the counter-phase flapping regime observed for flexible bodies placed at moderate distances [23, 24, 25].

The described dynamics is supported by Fig. 13, where the position of the three wings along with the pressure field at different instants within one flapping cycle is shown. From these snapshots, we argue that when two adjacent wings get closer (approximately when the pitching angle is at its maximum/minimum), the flow velocity has to increase due to the narrower effective cross section; hence, the pressure minimum gets amplified compared to the single case, this in turns increasing the amplitude of the lift force and, consequently, the wing oscillation.

In Fig. 12, the PP oscillation amplitude of W1 is compared to that of the isolated case. As mentioned before, the variation of this quantity does not provide a direct indication of the gain in performance that is achieved employing this configuration. To perform a more detailed comparison, let us focus on the quantities that are directly associated with the plunge power, whose time evolutions are reported in Fig. 14. The first one is the PP transverse velocity (Fig. 14a), for which the observed difference, in terms of amplitude, is about $20\%$ for all wings. The second is the lift coefficient (defined as $C_{L}=2\,F_{\mathrm{aero}}^{y}/(\rho_{\mathrm{f}}csU^{2})$), shown in Fig. 14b, for which the variation is approximately $20\%$ for W1 and $40\%$ for the side wings. For the latter, the increment is even more pronounced since the flapping motion is asymmetric, so that higher positive (or negative) peaks of $C_{L}$ are found for the top (or bottom) wing. By multiplying these two quantities, we obtain the instantaneous power coefficient reported in Fig. 14c, where one can see how the positive peaks in the side-by-side configuration are appreciably higher than in the single case, while the negative peaks remain essentially at the same values. Consequently, the average power coefficient undergoes the significant increase already presented in Fig. 11c.

In light of these results, we move further by considering arrays of more than three objects, fixing the mutual distance between adjacent devices to $r_{y}=1$, for which we found the highest increase in performance. The outcome of this analysis is presented by Fig. 15 collecting the average power coefficient for each device: when placing two additional plates at $y=\pm 2$ (i.e., considering an array of $N_{\mathrm{d}}=5$ devices), these achieve a further enhanced performance compared to the case where $N_{\mathrm{d}}=3$. The same trend holds when increasing the number of devices to $N_{\mathrm{d}}=11$. Increasing $N_{\mathrm{d}}$, the motion of flapping objects may be found to be less regular and periodic, as e.g. for the side wing placed at $y=5$. Hence, we now evaluate the average in a statistical sense, assessing its convergence in time (this corresponding to error bars in Fig. 15).

Emphasizing the comparison with the total power expected considering $N_{\mathrm{d}}$ isolated devices, we can write the overall power coefficient associated with the whole array $\bar{C}_{P}^{\mathrm{tot}}$ as:

[TABLE]

where $\bar{C}_{P}^{\mathrm{S}}$ is the power coefficient of the single device in isolated configuration and the interference quantity $I=I(N_{\mathrm{d}})$ has been introduced. The latter quantifies the additional power due to the cooperative effect. Similarly, the relative increment can be expressed by introducing the interference factor

[TABLE]

Table 1 collects the values obtained in our computational study, from which $\mathit{IF}$ seems to converge with the number of devices approximately to unity, i.e. the amount of available power is almost $100\%$ increased compared to the sum of ${N_{\mathrm{d}}}$ isolated devices. Moreover, this evidence can be representative when extrapolating our results to the situation of a network made of a vast number of devices. In this situation the side effects, although beneficial in terms of performance, are expected to be negligible. Nevertheless, each device can be thought to behave as those in the central part of the array. We can therefore conclude that employing a side-by-side arrangement of many devices appears to be substantially beneficial in terms of EH potential.

4.4.1 Wind-tunnel experiments

In order to corroborate the numerical findings, our study is complemented by an experimental analysis of the side-by-side arrangement made of an array of three devices, shown in Fig. 16. Experiments were performed using the wind-tunnel facility at Physics Department of the University of Genoa, whose characteristic parameters, along with the procedure used to acquire the motion of flapping wings in time, have been described in previous works [26, 28, 30].

Each of the three components is similar to the devices already presented in Refs. [28, 29] and is briefly described as follows. A $0.1\mathrm{mm}$-thick foil made of polyvinil acetate, with planform dimensions $c=30\mathrm{mm}$ and $s=65\mathrm{mm}$, is glued to a brass rod that is connected (allowing free rotation) to four elastomeric elements made of siliconic rubber, aligned with the streamwise direction and pre-stretched so that, following Ref. [28], the value of an effective stiffness can be estimated to be $\mathcal{K}^{\mathrm{eff}}_{(y)}\approx 13\mathrm{N/m}$. The mass of the moving body is $m=0.94\mathrm{g}$ and the center of mass is approximately placed $0.2c$ aft of the leading edge.

The three devices are collocated within a fixed frame and placed approximately at the center of the test chamber. The mutual distance between devices corresponds to $r_{y}\approx 1.1$. A freestream air velocity $U=4.5\mathrm{m/s}$ is chosen, at which regular LCO was found for devices tested in the single configuration.

Assuming air at standard conditions, the chord-based Reynolds number can be estimated as $\mathit{Re}\approx 9\times 10^{3}$, while from the other quantities we can derive the remaining equivalent nondimensional parameters: $\AR\approx 2$, $\rho_{\mathrm{w}}\approx 17$ and $K\approx 21$. Except for the aspect ratio, these values are clearly different from those considered in the numerical investigation. The analysis is thus intended as complementary, in order to assess the robustness of the outlined mechanism when moving into the operational range of the real EH application.

The experiment has been conducted as follows. First, we performed measurements of each device taken individually, i.e. removing the other two from the array and retaining the device under consideration in the same position as in the multiple configuration. Then, we tested the multiple configuration, where all three devices are present.

The time traces of the acquired PP oscillation for each device are reported in Fig. 17, where it can be observed that the behaviour is analogous to that of numerical simulations: the amplitude of the oscillation increases and motions of adjacent wings are essentially in counterphase, as it is also shown by the side views of Fig. 18 (see also Movie 3 in the Supplementary Material).

To quantify the effect due to the mutual interaction, Table 2 reports values of the same observables analysed numerically, i.e. amplitude of PP oscillation, Strouhal number and power coefficient, both for the single and multiple configurations. Focusing on data of the single configuration, we observe that a certain difference exists between the three devices, which can be ascribed mainly to constructive details. Nevertheless, when considering the same quantities in the multiple configuration, an increase with respect to the single case is found for all devices. In particular, the PP amplitude is maximised for the central wing (W1), with an increase around $40\%$ with respect to the individual configuration. Variation in frequency is smaller, i.e. between $2$ and $10\%$ with the same resulting synchronisation between devices that is observed numerically. From the acquired movies it is also possible to evaluate the aerodynamic forces and thus the power associated to the plunge and pitch motions. In Table 2, we report the cycle-averaged plunge power coefficient, as defined in Eq. (8). The increment when the three wings are coupled is evident, confirming qualitatively the results obtained by the simulations.

From the quantitative viewpoint, some differences can be noted between the experimental and the numerical results (i.e., comparing Table 2 with Fig. 3 for the single device and Fig. 11 for the side-by-side arrangement). We underline, however, that the experimental configuration and the idealized model differ in several features that may explain the observed discrepancies: for example, for the former we have: (i) higher $\mathit{Re}$ number, (ii) nonhomogeneous mass distribution and (iii) nonlinear elastomeric elements in place of Hooke springs. These features appear to influence the device dynamics and on-going investigations are aimed at characterizing more exhaustively their role. In particular, we plan to perform simulations for a configuration where such details are included, hence further more representative of the real EH device.

5 Conclusion

This work has investigated the dynamics of an aeroelastic system suitable for EH purposes, based on fluttering oscillations of elastically-bounded plates in a laminar flow – focusing, in particular, on the interaction between multiple devices arranged into three basic configurations (i.e.: in-line, staggered and side-by-side). As a first step, we have considered an effective physical model where homogeneous plates are anchored by linear springs and immersed in uniform flow. By employing a suitable numerical procedure based on a finite-difference Navier-Stokes solver coupled with a moving-least-squares immersed boundary method, we have performed numerical simulations for the single and multiple wing configurations. Focusing on the main observables of flapping motion such as the amplitude and frequency of oscillation, as well as performance-related quantities (i.e. power coefficient and Betz efficiency), we investigated how the dynamical behaviour is affected by the governing parameters (e.g., the nondimensional stiffness) and the mutual distance at which devices are placed in the multiple configuration.

For both the in-line and staggered arrangements, the performance of interacting devices is found to be worse than in the single case, although a recovery in performance can be achieved by tuning the elasticity of downstream devices. When considering the side-by-side configuration, in contrast, the interaction turns out to be beneficial and relevant increases of all quantities of interest are found. Further enhancements are obtained when increasing the number of devices in the array. A constructive interference is found, causing the total power coefficient of the network to increase up to $100\%$ with respect to the expected amount by the same number of isolated devices.

The occurrence of this cooperative effect has been verified experimentally in a case representative of the real EH application, considering an array of prototypal devices in a higher Reynolds flow, revealing good agreement with numerical findings despite the different configurations. The present results therefore suggest the development of dense arrays of flutter-based EH devices as a strategy that could enable higher performance for this kind of novel technology.

Acknowledgements

The authors acknowledge CINECA and INFN for the availability of high performance computing resources and support.

Appendix A Validation of numerical method

The numerical procedure used in the present study has been extensively validated for several FSI problems involving both rigid and deformable bodies [33, 35]. In the current investigation, however, we assess the dependency of numerical results on the spatial and temporal resolution. To this end, the chord-based Reynolds number is set to $\mathit{Re}=1000$, i.e. the highest value tested in the whole study, and we choose the baseline values for the other control parameters, i.e. $\AR=2$, $\rho_{\mathrm{w}}=2$ and $K=3$, this configuration corresponding to a regular flapping state of finite amplitude. A finite initial perturbation on the angular velocity is given in order to shorten the transient, although the achievement of the same steady limit-cycle was verified in case of unperturbed initial condition.

The wing is initially placed with its geometrical center at the origin $(0,0,0)$. The domain dimensions are set at $\pm 5$ both in the $x$ and $y$ direction, and at $-5$ and $10$ in the $z$ direction. A Cartesian grid is used, with uniform resolution in the spanwise direction and stretched in the other two directions. However, the grid spacing $h$ is made locally uniform in a region close to the wing, i.e. for $-1<z<2$ and $-1<y<1$. In order to quantify convergence with respect to the spatial resolution, eleven different grids were considered, characterized as shown in Table 3. While the same topology is retained, the grids differ by an overall refinement factor. The wing is assumed to have zero thickness and is thus discretized in a surface mesh made of triangular elements whose characteristic length is proportionally adjusted to be essentially equal to the minimum Eulerian spacing $h$.

The approximate solution can be expressed as a function of grid spacing as follows:

[TABLE]

where $h$ is the grid spacing, $f^{*}$ is the exact solution (i.e., the solution that one would have for $h\to 0$), $C$ is a constant and $p$ is the (actual) order of convergence. Instead of fitting the obtained data using this expression, we find more convenient to introduce the relative error with respect to the finest grid (for which $h=h_{\mathrm{min}}$), defined as

[TABLE]

since by combining Eqs. (12) and (13), the following relation can be found:

[TABLE]

where only $C$ and $p$ appear. The fitting is therefore applied to this latter expression, from which we later extrapolate the exact value $f^{*}$.

Grid convergence is assessed by considering the amplitude of the transverse oscillation of the pivot point $\mathrm{Amp}(y_{\mathit{PP}})$. Each case is computed up to measure at least 5 cycles of steady LCO. Fig. 19a shows the corresponding mean values along with the minimum and maximum ones. From the plot one can notice that the convergence trend has a change at about $h=0.02$, the actual order improving while increasing the resolution. By applying Eq. (14), the curve depicted in Fig. 19a is found, where the coefficients appearing in the expression are $C\approx-40.7$ and $p\approx 1.44$, yielding the extrapolated value of the exact solution $f^{*}\approx 0.56$. In the plot we also report data from two cases employing uniform grids with constant spacing $h=\{0.02,0.04\}$, showing that the discrepancy with equivalent stretched grids (with same $h$) looks contained and decreases while increasing the resolution. Moreover, we can also refer directly to the relative error with respect to the finest grid considered, as defined by Eq. (13). This quantity is shown in Fig. 19b where we can notice overall second-order accuracy for finer grids and first-order accuracy for coarser ones (in agreement with previous studies [33]).

In the presented parametric studies, two different grid settings are used. For the single (Sec. 4.1) and in-line (Sec. 4.2) configurations, grid E is chosen. For the staggered (Sec. 4.3) and side-by-side (Sec. 4.4) configurations, we use instead a uniform grid with $h=0.04$. The latter choice is motivated by the fact that a wider region of the domain has to be refined in this case. One can note that the oscillation gets underestimated when using coarser grids, thus the numerical result is generally conservative. A similar trend can be found for the pitching amplitude, while the flapping frequency shows a convergence that is far more rapid with variations less than $2\%$ within the considered range of spatial resolution.

Fixing the choice of the grid, we are able to investigate the dependence of the solution with respect to the simulation timestep. Fig. 20 shows the convergence when decreasing the timestep, with an estimated error of about 3% with respect to the extrapolated value for vanishing $\Delta t$. Last, the sensitivity to the domain size was checked, with negligible differences observed in the resulting flapping motion when doubling the size of the bounding box in all directions.

We conclude by providing a comparison of our results with those obtained with an essentially alternative approach, i.e. using a body-conforming mesh. The computation is performed using the open source library OpenFOAM [37], which is based on the finite volume method and offers dynamic mesh features, such as the overset treatment that is employed in our benchmark. For this test, we set the case so that the spatial and temporal resolutions are comparable to those of grid A. In Fig. 21, we report the time histories of the pivot point motion and pitching angle obtained with the two approaches: all quantities look in good agreement, with small differences in the amplitude and period of the oscillation.

The simulations were performed on the same workstation using 8 processors. The comparison between the computational times of the two codes indicates a wall clock time of $\simeq 0.25\mathrm{s}$ per time step for the IBM code and $\simeq 5\mathrm{s}$ for OpenFOAM, thus yielding a factor 20. Although the specific figures might depend on the particular computer architecture and the details of the problem, it is clear that the IBM code performs more efficiently for this class of problems.

It is worth mentioning that the differences might become even more evident when more than one flapping element is considered in the problem. In fact, while for IB methods the CPU time would grow only because of the larger number of immersed surface elements, for a moving grid method there would be more grid patches in relative motion and the computational load would increase more rapidly than the surface element counting.