Constructive Aerodynamic Interference in a Network of Weakly Coupled Flutter-Based Energy Harvesters

: Converting ﬂow-induced vibrations into electricity for low-power generation has received growing attention over the past few years. Aeroelastic phenomena, good candidates to yield high energy performance in renewable wind energy harvesting (EH) systems, can play a pivotal role in providing su ﬃ cient power for extended operation with little or no battery replacement. In this paper, a numerical model and a co-simulation approach have been developed to study a new EH device for power generation. We investigate the problem focusing on a weakly aerodynamically coupled ﬂutter-based EH system. It consists of two ﬂexible wings anchored by cantilevered beams with attached piezoelectric layers, undergoing nonlinear coupled bending–torsion limit cycle oscillations. Besides the development of individual EH devices, further issues are posed when considering multiple objects for realizing a network of devices and magnifying the extracted power due to nonlinear synergies and constructive interferences. This work investigates the e ﬀ ect of various external conditions and physical parameters on the performance of the piezoaeroelastic array of devices. From the viewpoint of applications, we are most concerned about whether an EH can generate su ﬃ cient power under a variable excitation. The results of this study can be used for the design and integration of low-energy wind generation technologies into buildings, bridges, and built-in sensor networks in aircraft structures.


Introduction
The energy harvesting concept is related to the small-scale energy extraction of various sources, usually by means of unconventional mechanisms. For example, the unavoidable random mechanical vibrations of a structure, such as a building, a bridge, or an aircraft wing, may be transformed into an electrical signal with a properly designed transducer. Depending on the application, this electrical signal could potentially replace small batteries to power small devices, such as wireless sensors and transmitters for health monitoring, microelectromechanical (MEM) devices, low maintenance, and hard-to-access actuators, to mention a few examples.
Rather than feeding on the vibrations of existing systems, a popular option for wind energy harvesting relies on external devices, harvesters, especially designed to exploit one or more fluid-interaction

Materials and Methods
In the following sections, the different components of the computational framework are presented. First, the aerodynamic model and the Unsteady Vortex Lattice Method is succinctly described. Next, the structural model and finite-element discretization are discussed. Finally, the coupling technique between these physics is laid out.

Aerodynamic Model
The main purpose of the aerodynamic model in the present framework is to provide an estimate of the aerodynamic loads. The well-known Unsteady Vortex Lattice Method (UVLM) was chosen for this purpose because it constitutes a trade-off between desired precision and computational cost compared to other flow solvers based on computational fluid dynamics (CFD) techniques such as Reynolds Average Navier-Stokes (RANS), Large Eddy Simulation (LES), Direct Numerical Simulation (DNS), and so on.
The UVLM is based on the premises of incompressible and irrotational flow outside the compact vorticity containing regions: boundary layers adjacent to the solid boundaries and wakes. In real flows at moderate to high Reynolds numbers, boundary layers develop adjacent to the surfaces of Aerospace 2020, 7, 167 4 of 29 the bodies. Enough vorticity is present in this layer to smoothly build up the flow velocity from zero at the surface (no-slipping condition) to the outer flow condition. As the Reynolds number increases, the thickness of the boundary layer decreases due to the progressive reduction of the influence of viscous diffusion relative to convective linear momentum transport. In the limit of infinite Reynolds number, the boundary layers (and other regions containing vorticity) become extremely thin but still enclose enough vorticity to produce the velocity transition from no-slip to outer flow conditions. In fact, in this situation, the boundary layers appear to produce an almost discontinuous jump in the tangential velocity of the flow, resembling the ideal scenario of an infinitesimally thin vortex sheet. Two types of vortex sheets are recognized: bound vortex sheets modeling the boundary layers attached to the solid boundaries ∂B and free-vortex sheets representing the wakes ∂F. For problems involving N B bodies, the solid boundaries ∂B i are collectively represented by ∂B = ∪ N B i=1 ∂B i . Similarly, all N F vortex sheets are collectively represented by ∂F = ∪ N F i=1 ∂F i . The bound vortex sheets move with the bodies; as a result, a finite pressure jump develops across them. Conversely, the free-vortex sheets grow and deform freely to attain force-free positions. The two types of vortex sheets merge at the points of the bodies where flow separation occurs. In this effort, it is assumed that the immersed bodies (harvesters) are very thin; therefore, the vortex sheets on both sides of the body effectively collapse into one.
The flow outside the vortex sheets is incompressible and irrotational; hence, a velocity potential Ψ exists and satisfies the Laplace equation (continuity for incompressible, irrotational flows): The velocity field V(r; t) = ∇Ψ must satisfy at all times the no-penetration condition on the surface ∂B of the bodies: (V(r; t) − V S (r; t)) ·n(r; t) = 0 ∀ r ∈ ∂B, wheren(r; t) is the unit vector normal to the solid surface at the location r and V S (r; t) is the velocity of the solid surface at r. Additionally, it is required that the perturbations introduced on the flow by the bodies decay to zero far away from them. This regularity condition at infinity is satisfied automatically by the fundamental solutions (vortex, source/sink, dipole) of the Laplace equation. Due to the linearity of the Laplace equation, the fluid velocity anywhere in the flow field can be decomposed into the free-stream component V ∞ and the fluid velocity due to the vorticity present on the bound and free-vortex sheet V ∂B and V ∂F , respectively: V(r; t) = V ∞ + V ∂B (r; t) + V ∂F (r; t).
Aerospace 2020, 7, 167 5 of 29 Finally, the unsteady Bernoulli's equation can be used to obtain the pressure distribution on the body surface and the resultant aerodynamic loads. This equation is given by: The resultant aerodynamic loads, force, and moment about a pivot point O i acting on each body B i are computed as follows: In the UVLM formulation, the bound vortex sheets attached to each body are discretized with collections of finite-length vortex segments of constant circulation ∆Γ. In order to fulfill the spatial conservation of vorticity (Equation (6)), the vortex segments, in turn, must conform to an array of (ideally four-sided) closed vortex rings with constant circulation (Figure 1). Similarly, the wakes are discretized with vortex segments forming closed vortex rings of constant circulation that automatically satisfy the spatial conservation of vorticity.
Aerospace 2020, 7, x FOR PEER REVIEW 5 of 28 The resultant aerodynamic loads, force, and moment about a pivot point acting on each body are computed as follows: In the UVLM formulation, the bound vortex sheets attached to each body are discretized with collections of finite-length vortex segments of constant circulation ΔΓ. In order to fulfill the spatial conservation of vorticity (Equation (6)), the vortex segments, in turn, must conform to an array of (ideally four-sided) closed vortex rings with constant circulation (Figure 1). Similarly, the wakes are discretized with vortex segments forming closed vortex rings of constant circulation that automatically satisfy the spatial conservation of vorticity.
A consequence of the discretization is that the no-penetration condition on a body can only be satisfied at a finite number of collocation points on its solid surface. It can be shown [36] that the contribution to the normal component of the fluid velocity at the location ∈ of the collocation point i on the body due to all the bound vortex sheets can be approximated as follows: Figure 1. Schematic representation of the vortex lattice discretization of the bound vortex sheet on a solid surface. Finite-length vortex segments with circulations ΔΓ 1 , and ΔΓ 2 , are used to discretize the continuous vorticity distribution over the surface. The subindices 1 and 2 indicate two directions of discretization. When connected, these vortex segments form vortex rings whose circulation is conserved spatially.
where denotes the aerodynamic influence of the j-th vortex ring attached to the body over the i-th collocation point on the surface of body , is the circulation of the j-th vortex ring attached to the body , and is the number of bound vortex rings in . The influence 2 are used to discretize the continuous vorticity distribution over the surface. The subindices 1 and 2 indicate two directions of discretization. When connected, these vortex segments form vortex rings whose circulation is conserved spatially.
A consequence of the discretization is that the no-penetration condition on a body B k can only be satisfied at a finite number M k of collocation points on its solid surface. It can be shown [36] that the contribution to the normal component of the fluid velocity at the location r k i ∈ ∂B k of the collocation point i on the body B k due to all the bound vortex sheets can be approximated as follows: Aerospace 2020, 7, 167 6 of 29 where A km ij denotes the aerodynamic influence of the j-th vortex ring attached to the body B m over the i-th collocation point on the surface of body B k , G m j is the circulation of the j-th vortex ring attached to the body B m , and N m is the number of bound vortex rings in B m . The influence coefficients A km ij are obtained by summing up the contributions to the normal component of the fluid velocity at the collocation point i due to the vortex segments (with unitary circulation) that make up the vortex ring j. The velocity due to a vortex segment of length L and unit circulation is given by: where quantities L, r 1 , r 2 ,ê 1 , andê 2 are defined in Figure 2.
Aerospace 2020, 7, x FOR PEER REVIEW 6 of 28 coefficients are obtained by summing up the contributions to the normal component of the fluid velocity at the collocation point i due to the vortex segments (with unitary circulation) that make up the vortex ring j. The velocity due to a vortex segment of length L and unit circulation is given by: where quantities , 1 , 2 , ̂1, and ̂2 are defined in Figure 2. It is observed that ( ) blows up when 1 ∥ ; that is to say, when the point where the velocity is evaluated lies along the direction of the vortex segment. This situation can potentially lead to numerical instabilities when updating the position of the free-vortex segment in the wakes. To avoid such a situation, it is customary to implement a smoothing or cut-off technique. In the current framework, a simple cut-off strategy is used that consists of replacing the denominator in Equation (11) by ‖ × 1 ‖ 2 2 + ( ) 2 [37], which effectively adds a finite core to the vortex segment. The nopenetration condition at the collocation point i of body is readily expressed as follows: It is noted that the contribution of the wakes (free-vortex sheets ) was moved to the righthand side of Equation (12). Once generated, the circulation of the free-vortex rings is conserved in time, and the rings move with the local fluid velocity to attain force-free configurations. Therefore, not only is their circulation known but also their location; therefore, their contribution to the fluid velocity on the collocation points is known through the Biot-Savart equation (Equation (5)).
It can be shown that the following algebraic system of equations is obtained upon enforcing the no-penetration condition on all the collocation points on all the bodies: , and RHS is the normal component of the fluid velocity on the i-th control point of body due to all other agents (free current, wakes, body motion, and so on) but the bound vorticity (vortex rings).
In Equation (13), a quasi-steady aerodynamic coupling between the different bodies is captured It is observed that V L (r) blows up when r 1 L; that is to say, when the point where the velocity is evaluated lies along the direction of the vortex segment. This situation can potentially lead to numerical instabilities when updating the position of the free-vortex segment in the wakes. To avoid such a situation, it is customary to implement a smoothing or cut-off technique. In the current framework, a simple cut-off strategy is used that consists of replacing the denominator in Equation (11) by L × r 1 2 2 + (δL) 2 [37], which effectively adds a finite core to the vortex segment. The no-penetration condition at the collocation point i of body B k is readily expressed as follows: It is noted that the contribution of the wakes (free-vortex sheets ∂F) was moved to the right-hand side of Equation (12). Once generated, the circulation of the free-vortex rings is conserved in time, and the rings move with the local fluid velocity to attain force-free configurations. Therefore, not only is their circulation known but also their location; therefore, their contribution to the fluid velocity on the collocation points is known through the Biot-Savart equation (Equation (5)).
It can be shown that the following algebraic system of equations is obtained upon enforcing the no-penetration condition on all the collocation points on all the bodies: , and RHS k i is the normal component of the fluid velocity on the i-th control point of body B k due to all other agents (free current, wakes, body motion, and so on) but the bound vorticity (vortex rings).
In Equation (13), a quasi-steady aerodynamic coupling between the different bodies is captured by the off-diagonal matrix blocks [A km ij ] with k m, which account for the bound vorticity interactions (aerodynamic shadow effect) of body B m over body B k . The other unsteady aerodynamic coupling between the bodies is encapsulated in the right-hand side terms RHS in which the effect of the wake of a body over another body is accounted for (Equation (12)).
The free-vortex sheets or wakes are created from the edges of the bodies where flow separation occurs by means of vorticity convection. Flow separation cannot be predicted within an irrotational framework (no-slip condition cannot be enforced) unless a boundary layer model or another more sophisticated strategy is incorporated. In the current effort, no such strategies are implemented. Instead, flow separation is assumed a priori to occur from specified locations, particularly from the trailing edge of streamlined bodies. The flow separation phenomenon is replicated in the current framework by the convection of the vortex segments along the trailing edge of the bodies into the flow field ( Figure 3). By Equation (6), the advected vortex segments must also form closed vortex rings, hence they close on the trailing edge through vortex segments parallel to the flow field. Once a part of the wake, the vortex rings move with the local fluid velocity conserving its sectional circulation. The motion of a free-vortex ring is achieved by moving its four corner points with the local fluid velocity. The time stepping is done through the first-order forward Euler scheme as follows: Aerospace 2020, 7, x FOR PEER REVIEW 7 of 28 aerodynamic coupling between the bodies is encapsulated in the right-hand side terms RHS in which the effect of the wake of a body over another body is accounted for (Equation (12)). The free-vortex sheets or wakes are created from the edges of the bodies where flow separation occurs by means of vorticity convection. Flow separation cannot be predicted within an irrotational framework (no-slip condition cannot be enforced) unless a boundary layer model or another more sophisticated strategy is incorporated. In the current effort, no such strategies are implemented. Instead, flow separation is assumed a priori to occur from specified locations, particularly from the trailing edge of streamlined bodies. The flow separation phenomenon is replicated in the current framework by the convection of the vortex segments along the trailing edge of the bodies into the flow field ( Figure 3). By Equation (6), the advected vortex segments must also form closed vortex rings, hence they close on the trailing edge through vortex segments parallel to the flow field. Once a part of the wake, the vortex rings move with the local fluid velocity conserving its sectional circulation. The motion of a free-vortex ring is achieved by moving its four corner points with the local fluid velocity. The time stepping is done through the first-order forward Euler scheme as follows:

Structural Model
The main aspects of the structural model are summarized next. The equations of the motion of a finite-element, three-dimensional piezoelectric beam and its associated set of reduced-order equations are presented.
In this effort, the piezoelectric beam studied is composed of an elastic core and two exterior

Structural Model
The main aspects of the structural model are summarized next. The equations of the motion of a finite-element, three-dimensional piezoelectric beam and its associated set of reduced-order equations are presented.
In this effort, the piezoelectric beam studied is composed of an elastic core and two exterior layers of piezoelectric material (transducers) with their respective electrodes. The piezoelectric layers are symmetrically arranged with respect to the bending neutral axis of the section. This beam is depicted schematically in Figure 4a.
Aerospace 2020, 7, x FOR PEER REVIEW 8 of 28 accounted for. Therefore, two generalized coordinates, and , are employed for the electric potential difference across the upper and lower piezoelectric layer, respectively. To estimate the power, the electrodes in each layer are connected to the resistors and , as shown in Figure 4b. This electric configuration has been extensively used in similar studies.
Regarding the material properties, the elastic layer is assumed to be homogeneous and isotropic. On the other hand, the piezoelectric layers are idealized as transversely isotropic materials, with their isotropic plane normal to the poling direction [38]. In this effort, a Rayleigh-like beam element satisfying the following kinematic hypothesis is assumed: H1. The cross-section remains planar and does not suffer in-plane deformations. This means that in-plane shear strains and Poisson's effect due to axial strains are assumed negligible; H2. The rotated cross-sections are always orthogonal to the deformed beam axis. This is the well-known Euler-Bernoulli assumption; and H3. The elastic core and the piezoelectric layers are perfectly bonded. This means that the deformation field is continuous across the elastic and piezoelectric layers.
In the description that follows, Latin indices i, j, k, and so on range from 1 to 3, and Greek indices range from 2 to 3. In addition, Einstein's summation is implied unless otherwise specified.
The body kinematics under deformation are described by means of a reference and the current configurations, ℬ and ℬ , respectively, both defined with respect to a global reference frame = {̂1,̂2,̂3}. In order to describe the reference configuration, a right-handed orthonormal frame ℰ = {̂1,̂2,̂3}, called reference frame, is introduced, with an associated set of coordinates ( 1 , 2 , 3 ). This frame is such that ̂2 and ̂3 lay on a generic cross-section and ̂1 is parallel to the beam axis. The beam reference configuration is described by the position vector 0 =̂. To describe the beam current configuration, another right-handed orthonormal frame ℯ = {̂1,̂2,̂3}, called moving or current frame, is introduced, such that ̂= ( 1 )̂ and ( 1 ) ∈ ( ) is a one-parameter rotation tensor. Due to the hypothesis (H1), the position in the deformed configuration of material points lying on an arbitrary cross-section are uniquely specified by the vector ( 1 , 2 , 3 ( 1 ) is the position in the current configuration of a point lying on the beam axis and =̂. The beam reference domain B is decomposed as the Cartesian product of the beam axis P and the beam normal cross-section A. The cross-section A is decomposed in three subdomains A i ; that is to say, where A 1 and A 3 are the cross-sections of the piezoelectric layers and A 2 is the cross-section of the elastic layer. Each cross-section A i is rectangular with the same width b and thickness t i (Figure 4a). Each piezoelectric layer is covered by a pair of electrodes on their top and bottom surfaces. These electrodes are assumed to be perfectly conductive and to have negligible mechanical properties ( Figure 4c). In addition, the transducers are electrically independent of one another, which means that the coupling between the axial and bending deformations can be accounted for. Therefore, two generalized coordinates, v u e and v l e , are employed for the electric potential difference across the upper and lower piezoelectric layer, respectively. To estimate the power, the electrodes in each layer are connected to the resistors R u and R l , as shown in Figure 4b. This electric configuration has been extensively used in similar studies.
Regarding the material properties, the elastic layer is assumed to be homogeneous and isotropic. On the other hand, the piezoelectric layers are idealized as transversely isotropic materials, with their isotropic plane normal to the poling direction [38].
In this effort, a Rayleigh-like beam element satisfying the following kinematic hypothesis is assumed: The cross-section A remains planar and does not suffer in-plane deformations. This means that in-plane shear strains and Poisson's effect due to axial strains are assumed negligible; H2. The rotated cross-sections are always orthogonal to the deformed beam axis. This is the well-known Euler-Bernoulli assumption; and Aerospace 2020, 7, 167 9 of 29 H3. The elastic core and the piezoelectric layers are perfectly bonded. This means that the deformation field is continuous across the elastic and piezoelectric layers.
In the description that follows, Latin indices i, j, k, and so on range from 1 to 3, and Greek indices range from 2 to 3. In addition, Einstein's summation is implied unless otherwise specified.
The body kinematics under deformation are described by means of a reference and the current configurations, B and B D , respectively, both defined with respect to a global reference frame N = n 1 ,n 2 ,n 3 . In order to describe the reference configuration, a right-handed orthonormal frame E = Ê 1 ,Ê 2 ,Ê 3 , called reference frame, is introduced, with an associated set of coordinates X 1 , X 2 , X 3 . This frame is such thatÊ 2 andÊ 3 lay on a generic cross-section andÊ 1 is parallel to the beam axis. The beam reference configuration is described by the position vector R 0 = X iÊ i . To describe the beam current configuration, another right-handed orthonormal frame e = ê 1 ,ê 2 ,ê 3 , called moving or current frame, is introduced, such thatê k = Λ X 1 Ê k and Λ X 1 ∈ SO(3) is a one-parameter rotation tensor. Due to the hypothesis (H1), the position in the deformed configuration of material points lying on an arbitrary cross-section are uniquely specified by the vector where r X 1 is the position in the current configuration of a point lying on the beam axis and X = X αÊ α .
In the present effort, only linear beam theory is considered. This means that the components of the displacement vector of a material point on the beam axis u 0 = r X 1 − X 1Ê 1 are small compared to the beam length L, and the rotations are small enough so that the rotation tensor can be approximated to the first order by Λ X 1 ≈ 1 +θ X 1 [40]. Here, 1 stands for the second-order identity tensor, andθ is a skew-symmetric tensor whose components can be interpreted as infinitesimal rotations. It can be shown that two successive infinitesimal rotations are independent of the order (commutative property), and thus, they can be regarded as vectors. With these assumptions, the displacement vector field of a generic point belonging to the beam cross-section can be expressed as, Here, θ 2 (θ 3 ) denotes the rotation of the cross-section around theÊ 2 (Ê 3 ) axis, and the axial vector associated toθ (2,3) is given by θ (2,3) = θ αÊ α . It is recalled at this time that the theory above only relates to the bending motion of the beam, so the rotation about the axisÊ 1 (torsion) is not included in order to fulfil Euler-Bernoulli hypothesis (H2) [41]. This is equivalent to the assumption that the torsional and bending deformation of the beam are independent (uncoupled). As it is customary in linear 3D beam formulations, the torsion effect is included in an ad-hoc fashion during the finite-element discretization procedure. Thus, a three-dimensional Euler-Bernoulli beam element is obtained, which includes torsion in a completely decoupled way.
From Equation (15), knowing that R = u + R 0 , the deformation gradient tensor F = ∂ k (R) ⊗Ê k can be computed, where the symbol ∂ k (·) stands for partial derivative with respect to X k and ⊗ represents the tensor product. The Green-Lagrange strain tensor ε GL can then be computed by means of the displacement gradient tensorF = F − 1. Neglecting all higher-order strain terms in ε GL , except Aerospace 2020, 7, 167 10 of 29 for the so-called von Kármán nonlinear strains, which involve the squares of ∂ 1 u 2 0 and ∂ 1 u 3 0 , leads to the following expression of the Green-Lagrange strain tensor [42], Due to the hypothesis (H2), the shear strain components (coordinates ofÊ 1 ⊗Ê 2 andÊ 1 ⊗Ê 3 ) in Equation (17) must be identically zero. This hypothesis leads to two additional kinematical constraints: 3) , or equivalently, κ = ∂ 1 θ αÊ α , can be recognized as the curvature vector of the beam axis.
The stress-strain constitutive relations for the elastic and piezoelectric layers are assumed to be linearly elastic, i.e., they follow Hooke's law. In addition, piezoelectric materials generate an electric field when they are deformed; conversely, they suffer deformation when an electric field is applied across them. These phenomena are called direct and converse piezoelectric effects, respectively [38]. The piezoelectric behavior can be described by linear constitutive equations if the electric fields and the mechanical strains are sufficiently small, yielding the following standard form where σ is the mechanical stress tensor, D is the electric displacement vector (or electric flux density), E e is the electric vector field, which is related to the electric potential v e as E e = −∇v e , C E e is the four-order elasticity tensor evaluated at a constant electric field, m is the third-order tensor of piezoelectric coefficients, e ε m is the second-order dielectric permittivity tensor evaluated at constant strain field, and the notation ":" indicates tensor double-contraction operation. It should be noted that the constitutive relations exposed in Equation (18) reduce to the well-known generalized Hooke's law σ = C : ε GL for elastic layers when the piezoelectric effect is ignored (E e ≡ 0).
As can be observed in Figure 4c, the piezoelectric layer is poled across its thickness. In addition, the electric potential is assumed to vary linearly across the electrode pair, namely, the electric field is uniform along the thickness direction. Therefore, the electric vector field through the i-th piezoelectric layer has only one nonzero component and can be expressed as the voltage-to-thickness To obtain the equations of motion for the electromechanical beam introduced above, an extended version of the variational Hamilton's principle that accounts for piezoelectric media is used [43,44]. There are two alternatives for dealing with the independent electrical variables. The first option is to use the electric charges as independent variables. The second option is to use the flux linkages as the independent variables. In any case, the electric fluxes can be related to the voltage through the constitutive relations of each electrical component. The resulting variational principle may be stated in Mandel's notation as follows [45], where the over-dots denote time derivatives, ρ is the density of the piezoelectric/elastic layer,ε ∈ R 6 is the Mandel representation of the second-order tensor ε GL , C E e ∈ R 6×6 is the Mandel representation of C E e , d m ∈ R 3×6 is the Mandel representation of m , and δW nc is the nonconservative virtual work, which in the present effort is done by the fluid forces (Equations (8) and (9)) and by the electric dissipation "forces." Due to the beam hypotheses adopted, ε 1 is the only nonzero component of the vectorε and therefore only C E e 11 is needed. Because of the electric field direction and the assumption of transversely isotropic material for the piezoelectric layers, matrices d m and e ε m are reduced to the scalars (d m ) 31 and (e ε m ) 33 , respectively.

Discretized Dynamic Formulation (Finite-Element Model)
The resulting partial differential equations obtained through the generalized Hamilton's principle are then discretized in the spatial domain by the finite-element method. This technique requires the beam domain to be divided into a finite number, N e , of subdomains called finite elements. In this work, a two-node first-order finite-element is adopted. The local element reference frame and associated Cartesian coordinates system ξ 1 , ξ 2 , ξ 3 are shown in Figure 5. The origin is located at the first node andÊ k 1 is parallel to the beam element axis, which is known as the elastic axis. Next, the displacement field given in Equation (15) is approximated by using a set of interpolation functions. In particular, Hermite cubic interpolation functions, H j (j = 1, 2, 3, 4), are adopted for the transversal displacement, and linear Lagrange interpolation functions, L j ( j = 1, 2), are used for the axial displacement. At this stage, the torsion effect is incorporated in an ad-hoc fashion and is also interpolated by using the linear functions L j [46]. Consequently, the displacement field u X 1 , X 2 , X 3 ; t is approximated, at the element level, as u k ξ 1 , ξ 2 , ξ 3 ; t = H k q k e . Here, H k ξ 1 , ξ 2 , ξ 3 is the so-called shape function matrix, and q k e = q k 1 , q k 2 T is the vector of nodal degree of freedom for the k-th element, where, Aerospace 2020, 7, x FOR PEER REVIEW 11 of 28 ∈ ℝ are the linear and nonlinear part of the electromechanical coupling vector, ( ) is the local nodal force vector, and and are the capacitance and the resistance of the s-th electric circuit ( = for top circuit and = for bottom circuit). For a purely elastic beam element, without piezoelectric layers, Equation (21) reduces to ̈+ [ + ( )] = ( ). The element geometric stiffness matrix and the nonlinear coupling terms arise from the von Kármán nonlinear terms included in the strain Green-Lagrange tensor along with the following two assumptions: (1) the (̂2 ,̂3 ) directions are principal axes of inertia, and (2) the axial forces are constant within each element. It is important to note that these axial forces and the nonlinear coupling terms must be calculated at every time step of the time-stepping procedure because they both depend on the configuration.
Finally, the global governing equations for a multifunctional structure made up of purely elastic and piezoelectric members are obtained through a standard FEM assembling procedure as follows where , , , ̅ and ̅ are the global mass matrix, elastic stiffness matrix, geometric stiffness matrix, and linear and nonlinear coupling matrices, respectively,  An adequate representation of an energy harvesting structure may require a large number of beam elements, which, in turn, leads to a large number of both nodal variables and degrees of freedom. Next, the modal decomposition method is used to reduce the dimension of the problem. This method can be performed in the following steps: (i) Obtain the structure's natural frequencies and eigenvectors (mode shapes): a. Express the nodal degrees of freedom as = . b. Replace this expression into the linear, undamped, and homogeneous version of Equation (22) to obtain the generalized eigenvalue problem ( − 2 ) = . c. Solve the eigenvalue problem by obtaining the natural frequencies and associated eigenvectors (free-vibration modes).
(ii) Obtain the actual response: After some algebraic manipulations, the following semidiscrete set of nonlinear ordinary differential equations (ODEs) for the k-th piezoelectric beam element in the structural mesh is obtained, where M e ∈ R 12×12 is the local mass matrix, K el e ∈ R 12×12 is the local elastic stiffness matrix, K geo e ∈ R 12×12 is the local geometric stiffness matrix, which depends on the configuration, Θ s Lin ∈ R 12 and Θ s NLin ∈ R 12 are the linear and nonlinear part of the electromechanical coupling vector, F e (t) is the local nodal force vector, and C s and R s are the capacitance and the resistance of the s-th electric circuit (s = u for top circuit and s = l for bottom circuit). For a purely elastic beam element, without piezoelectric layers, Equation (21)  q k e q k e = F e (t). The element geometric stiffness matrix K geo e and the nonlinear coupling terms Θ s NLin arise from the von Kármán nonlinear terms included in the strain Green-Lagrange tensor along with the following two assumptions: (1) the  directions are principal axes of inertia, and (2) the axial forces are constant within each element. It is important to note that these axial forces and the nonlinear coupling terms must be calculated at every time step of the time-stepping procedure because they both depend on the configuration.
Finally, the global governing equations for a multifunctional structure made up of purely elastic and piezoelectric members are obtained through a standard FEM assembling procedure as follows where M, K el , K geo ,Θ s Lin andΘ s NLin are the global mass matrix, elastic stiffness matrix, geometric stiffness matrix, and linear and nonlinear coupling matrices, respectively, is the global vector of external forces due to the surrounding fluid flow, v u e (v l e ) is the vector of generalized voltages, and q is the vector of generalized coordinates. Structural damping is introduced in the equations of motion through the Rayleigh proportional damping C d = αM + βK el , with α and β constant coefficients. Further details of the structural model can be found in [47].
An adequate representation of an energy harvesting structure may require a large number of beam elements, which, in turn, leads to a large number of both nodal variables and degrees of freedom. Next, the modal decomposition method is used to reduce the dimension of the problem. This method can be performed in the following steps: (i) Obtain the structure's natural frequencies and eigenvectors (mode shapes): a.
Express the nodal degrees of freedom as q = φ i e jω i t . b.
Replace this expression into the linear, undamped, and homogeneous version of Equation (22) to obtain the generalized eigenvalue problem Solve the eigenvalue problem by obtaining the natural frequencies ω i and associated eigenvectors φ i (free-vibration modes).
(ii) Obtain the actual response: a.
Express the nodal degrees of freedom as a linear combination of the free-vibration modes; that is to say, q = Φη, where η ∈ R N nm ×1 is a column vector whose components are the modal coordinates, and Φ = [φ 1 φ 2 |· · · |φ N nm ] ∈ R (12×N e )×N nm is the modal matrix, whose columns are the eigenvectors φ i such that Replace this expression into Equation (22) and solve for the modal weights η. c.
Recover the vector of physical coordinates q from the modal coordinates η. In Step (ii) b., the resulting equations of motion are ..
where Λ ω = Φ T K el Φ = diag ω 2 i is a diagonal matrix whose N nm diagonal elements are the squares of the natural frequencies ω i , (23) is a set of 3 × N nm ordinary differential equations for the modal coordinates η. The modal coordinate transformation implies a significant reduction in the dimension of the problem, from (3 × 12 × N e ) degrees of freedom to (3 × N nm ).

Communication between Models and Numerical Integration
A fundamental piece in the aeroelastic simulations is the strategy for transferring information (TI) between the aerodynamic grid and structural mesh [48]. This two-way information exchange consists of transferring forces from the aerodynamic grid to the structural mesh and displacements and velocities from the structural mesh to the aerodynamic grid. It is possible to relate the displacement of arbitrary points in the aerodynamic grid with the generalized displacements of some nodes in the structural mesh by means of the following expression where u A (t) is a column vector with dimension 3N m , which contains the three components of the displacement of the selected points of the aerodynamic grid, and q(t) is a column vector with dimension 6(N e + 1), which contains the components of generalized displacement of a group of nodes of the structural grid. The linear-interpolation mapping G AS : q → u A depends on: (i) the geometry of both grids; (ii) the selected points in the aerodynamic grid; and (iii) the particular type of finite elements used to discretize the structure. It is recalled that the aerodynamic loads act on the control points of the aerodynamic grid, which do not necessarily coincide with the nodes of the structural mesh. According to Preidikman [48], G AS is a block-matrix whose individual blocks relate the translations of a point B on the aerodynamic grid to the generalized displacements of a set of nodal points on the structural mesh. To materialize this relation, point B is linked to the elastic-axis of the beam by means of a rigid arm which is required to be contained in the plane defined byÊ 2 andÊ 3 . In this sense, this procedure can be thought of as a master-slave scheme, with the structural nodes being the master nodes and the aerodynamic control points the slave nodes ( Figure 6). Considering point A on the elastic axis in-between the structural nodes I and J, the displacement vector of the point B can be obtained as u B = u A X 1 ; t +θ X 1 ; t X * , where X * is the position vector of point B relative to point A, andθ is a skew-symmetric tensor, which accounts for immaterial rotations around axesÊ k (already defined in Section 2.2). The displacement vector of point B can then be approximated, at the element level, by u B = H k q k e . Finally, it can be shown that matrix G AS is a block-matrix containing H k -like submatrices. body m, invoking the arbitrariness of virtual displacements, and performing some algebraic manipulations lead to the following force transfer map, where | is the transpose of the linear-interpolation mapping evaluated at the control points. Further details of the TI, as well as examples of matrices, can be found in [48]. In the approach followed here, the electromechanical FEM and airflow UVLM models are treated as two separate, though not independent, subsystems of a single dynamical system. Between these two subsystems, information is exchanged bidirectionally in an iterative sequence in order to continuously improve the estimation of the overall system response until convergence is reached. The time-marching numerical scheme used for the UVLM (Equation (14)) is well-known and has been previously reported in the literature [48][49][50]. On the other hand, the time-marching procedure adopted for FEM of the harvester array is based on Hamming's fourth-order predictor-corrector method [51]. This method further requires the set of second-order differential equations presented in Equation (23) to be recast as a system of first-order ODEs.

Co-Simulation Strategy
Although the fluid (air) flow and piezoelastic structure are two physical fields that are independently modeled and computationally implemented of one another, the coupling procedure is indeed considered strong in the sense that information is bidirectionally exchanged, and the chosen time step, which advances the solution in time, is unique and shared by both simulation environments. The interested reader can consult references [52] and [33] for further details about the co-simulation strategy employed in the present effort.
The co-simulation strategy to compute the solution at time + Δ is implemented through the following steps: 1. The configurations of the wakes are updated. The new position of each node of the vortex segments in the wakes is estimated within the UVLM module according to Equation (14). Throughout the rest of the procedure for the current time step, and until convergence is achieved, the wake configuration remains frozen. 2. The current fluid flow variables (velocity field, pressure distribution, and aerodynamic loads) are solved within the UVLM module. The current aerodynamic loads ( + Δ ) computed at To relate the forces F(t) (on the right-hand-side of Equation (22)) that act upon the structure with the aerodynamic forces F A (t) given in Equation (8), an "abstract" equivalence is established between them. Concretely, it is required that the virtual work done by both system of forces under an arbitrary virtual displacement is the same, δW A = δW S . Considering the N m control points of body m, invoking the arbitrariness of virtual displacements, and performing some algebraic manipulations lead to the following force transfer map, where G T AS CP is the transpose of the linear-interpolation mapping evaluated at the control points. Further details of the TI, as well as examples of G AS matrices, can be found in [48].
In the approach followed here, the electromechanical FEM and airflow UVLM models are treated as two separate, though not independent, subsystems of a single dynamical system. Between these two subsystems, information is exchanged bidirectionally in an iterative sequence in order to continuously improve the estimation of the overall system response until convergence is reached. The time-marching numerical scheme used for the UVLM (Equation (14)) is well-known and has been previously reported in the literature [48][49][50]. On the other hand, the time-marching procedure adopted for FEM of the harvester array is based on Hamming's fourth-order predictor-corrector method [51]. This method further requires the set of second-order differential equations presented in Equation (23) to be recast as a system of first-order ODEs.

Co-Simulation Strategy
Although the fluid (air) flow and piezoelastic structure are two physical fields that are independently modeled and computationally implemented of one another, the coupling procedure is indeed considered strong in the sense that information is bidirectionally exchanged, and the chosen time step, which advances the solution in time, is unique and shared by both simulation environments. The interested reader can consult references [52] and [33] for further details about the co-simulation strategy employed in the present effort.
The co-simulation strategy to compute the solution at time t + ∆t is implemented through the following steps: 1.
The configurations of the wakes are updated. The new position of each node of the vortex segments in the wakes is estimated within the UVLM module according to Equation (14). Throughout the rest of the procedure for the current time step, and until convergence is achieved, the wake configuration remains frozen.

2.
The current fluid flow variables (velocity field, pressure distribution, and aerodynamic loads) are solved within the UVLM module. The current aerodynamic loads F A (t + ∆t) computed at the control points of the aerodynamic grid are transferred to the nodal points of the structural FEM mesh by means of Equation (25). With F(t + ∆t), the response of each harvester in the array is computed through Equation (23).

3.
The computed state of the structure (nodal displacements and velocities) is transferred to the aerodynamic mesh through Equation (22). With the updated configuration a new prediction or correction of the aerodynamic loads is obtained. 4.
Steps 2 and 3 are repeated until a convergence criterion is met. This criterion requires that the infinite norm of the relative error of the computed displacements between two consecutive iterations is less than 10 −8 .

5.
With the converged configuration of the structure, Step 3 is repeated one last time to obtain the final (converged) estimate of the aerodynamic loads.
It has been observed through a large number of simulations that, on average, between 2 and 4 iterations per time step are required when one harvester is considered and between 4 and 5 iterations when two harvesters are considered.

Numerical Results
The results obtained with the computational implementation of the methods described before are summarized in this section. First, in Section 3.1, a brief "validation" of the code is presented based on the case of a single isolated piezoelectric harvester immersed in a uniform free flow. This case has been previously studied and reported in the literature.
In Section 3.2, a comprehensive study on a two-harvester array is performed. In this case, the focus is on the effect of the mutual aerodynamic interaction between the harvesters on the overall behavior of the system and its energy harvesting capabilities.

Verification of the Numerical Model
In order to verify that the computational implementation is correct and able to reproduce physical phenomena, the study conducted by De Marqui et al. [53,54] Table 1.   The main results of the present analysis are summarized in Table 2. The flutter speed = 41.505 m/s of the system was predicted, slightly overestimating the value of 40.5 m/s reported by De Marqui et al.
To estimate the flutter speed, the bisection approach is employed. Simulations at increasing values of the wind speed are performed, and the system response is analyzed to detect bifurcation in Figure 7. Reproduction of the physical model proposed in [53,54]. In the context of the present framework, the equivalent beam model is discretized with a total of 10 finite elements: three coupled electromechanical beam elements are used for the portion where the piezoelectric transducers are located (Section A) and seven elastic beam elements for the remainder of the structure (Section B). Although 60 different structural modes and natural frequencies are obtained through modal analysis, the accuracy of the high-frequency modes is known to be poor. Therefore, to obtain the reduced set of equations of motions (Equation (23)), only the first five structural modes are used for the modal decomposition. Regarding the aerodynamic model, a lattice of 30 spanwise by 8 chordwise elements is used (Figure 8). The imposed free flow direction is parallel to the y-axis (cross-flow), as shown schematically in Figure 7.
The main results of the present analysis are summarized in Table 2   To estimate the flutter speed, the bisection approach is employed. Simulations at increasing values of the wind speed are performed, and the system response is analyzed to detect bifurcation in  To estimate the flutter speed, the bisection approach is employed. Simulations at increasing values of the wind speed are performed, and the system response is analyzed to detect bifurcation in the system response from a stable focus (damped oscillations) to a finite-amplitude limit cycle oscillation. Once this change is identified, simulations with successively smaller wind velocity intervals are performed until close enough lower and upper bounds of the critical flutter speed are identified to the desired precision. In this article, the upper bound so identified is referred to as "flutter speed." To understand the behavior of the system when the airflow speed is equal to the flutter speed, a transient simulation was performed. In this case, an initial condition of 3 deg pitch at the wing tip was imposed to trip the flutter response. In Figure 9, the time history of the tip displacement u tip (Figure 9a) and rotation θ tip (Figure 9b) obtained at flutter speed are presented. From the figure, it is possible to appreciate the presence of an oscillatory motion of constant amplitude, which corresponds to a limit cycle oscillation LCO (the phase portrait inserts in Figure 9a,b) associated with the nonlinearity of the aerodynamic model. the system response from a stable focus (damped oscillations) to a finite-amplitude limit cycle oscillation. Once this change is identified, simulations with successively smaller wind velocity intervals are performed until close enough lower and upper bounds of the critical flutter speed are identified to the desired precision. In this article, the upper bound so identified is referred to as "flutter speed." To understand the behavior of the system when the airflow speed is equal to the flutter speed, a transient simulation was performed. In this case, an initial condition of 3 deg pitch at the wing tip was imposed to trip the flutter response. In Figure 9, the time history of the tip displacement ( Figure 9a) and rotation (Figure 9b) obtained at flutter speed are presented. From the figure, it is possible to appreciate the presence of an oscillatory motion of constant amplitude, which corresponds to a limit cycle oscillation LCO (the phase portrait inserts in Figure 9a,b) associated with the nonlinearity of the aerodynamic model. In order to determine the flutter frequency, the Fast Fourier Transform (FFT) is applied to the ( ) and ( ) signals. The normalized FFT amplitudes are shown in Figure 10, along with the first five structural natural frequencies (i.e., free vibration frequencies) of the system ( = 1, … ,5). In order to determine the flutter frequency, the Fast Fourier Transform (FFT) is applied to the u tip (t) and θ tip (t) signals. The normalized FFT amplitudes are shown in Figure 10, along with the first five structural natural frequencies (i.e., free vibration frequencies) of the system ω i (i = 1, . . . , 5). In order to determine the flutter frequency, the Fast Fourier Transform (FFT) is applied to the ( ) and ( ) signals. The normalized FFT amplitudes are shown in Figure 10, along with the first five structural natural frequencies (i.e., free vibration frequencies) of the system ( = 1, … ,5). From Figure 10, it can be appreciated that at flutter speed, all frequency components of the response coalesce to a unique value: the so-called the flutter frequency ω f . In this particular case, it is found that ω f = 72.6 rad/s, which lies between the second bending frequency (ω 2 ) and the first torsion frequency (ω 3 ). Additionally, the calculated value of ω f is in good agreement with the value reported by De Marqui (ω DM f = 72.07 rad/s). Table 2 summarizes the comparison between the results obtained with the present numerical model and those from De Marqui. The first five structural natural frequencies predicted through the present beam model are in good agreement with those obtained by De Marqui et al. using their plate model. The discrepancies between these predictions do not exceed 3%. Furthermore, the flutter frequency and speed are also well approximated through the approach proposed in this effort. These results show that the simulation framework is appropriate to study the aeroelastic behavior of the piezoelectric-based power-generating wing.

Case Study: Two Aerodynamically Coupled Harvesters
The case study considered herein consists of an array of two harvesters stacked vertically and separated by a vertical distance h (see Figure 11a). The main goal of this study is to characterize the combined piezoaeroelastic behavior of the system, quantifying the effect of the aerodynamic coupling between the harvesters. In order to assess the aerodynamic interaction, the flutter onset wind speed and the electrical output of the two-harvester system are compared with that of an isolated harvester. The model consists of two identical cantilevered harvesters, each with identical properties and discretization as the ones used in Section 3.1. Each harvester produces a vortex wake that is shed from the trailing edge and from the free tip of the plate, as shown schematically in Figure 11b. discretization as the ones used in Section 3.1. Each harvester produces a vortex wake that is shed from the trailing edge and from the free tip of the plate, as shown schematically in Figure 11b.
For the first study on this configuration (Section 3.2.1), the influence of the separation distance h on the mutual interaction among harvesters and the global system characteristics are analyzed. For the last study (Section 3.2.2), the system's postcritical behavior is studied for a particular fixed value of the separation h. For the first study on this configuration (Section 3.2.1), the influence of the separation distance h on the mutual interaction among harvesters and the global system characteristics are analyzed. For the last study (Section 3.2.2), the system's postcritical behavior is studied for a particular fixed value of the separation h.

Influence of the Separation Distance
The flutter speed of the harvester array V A F was computed for increasing values of the separation distance h. The variation of the nondimensional flutter speed V A F with the vertical distance between the harvesters is shown in Figure 12 Here,

Influence of the Separation Distance
The flutter speed of the harvester array was computed for increasing values of the separation distance ℎ . The variation of the nondimensional flutter speed ̅̅̅̅ with the vertical distance between the harvesters is shown in Figure 12 Here, ̅̅̅̅ = / is the normalized flutter speed of the array with respect to the flutter speed = 41.505 m/s of an isolated harvester. Based only on the dependence of the flutter speed with the spacing h, two zones, favorable and unfavorable, are identified in Figure 12. It is recalled that each individual harvester in the system only undergoes self-sustaining limit cycle oscillations at or above the flutter speed. Therefore, a favorable characterization is given to the region where lower wind speeds are needed to achieve flutter. A lower flutter critical speed of the system allows the harvesters to operate over a wider range of freestream wind velocities and to exploit the flutter phenomena at lower wind speeds, something otherwise impossible for the single isolated harvester configuration. Therefore, this scenario is deemed advantageous.
As it is shown in Figure 12, the transition between the favorable and unfavorable interaction zones occurs at a value of spacing ℎ⁄ approximately equal to five. In the favorable zone, the flutter speed of the array is below the flutter speed of an isolated harvester ( ̅̅̅̅ < 1). Conversely, in the unfavourable zone, the flutter speed of the array is above the flutter speed of an isolated harvester ( ̅̅̅̅ < 1). In this last zone, it is observed that a maximum value of ̅ is reached approximately at ℎ⁄ = 7.5. From this value on, on decreases asymptotically approaching the reference value ( ̅̅̅̅ = 1). This implies that as the distance increases beyond a certain point, the mutual aerodynamic interaction between harvesters becomes insignificant, and each harvester behaves effectively as if it was isolated. Continuing with the analysis, the amplitudes of the LCO of the wing tip pitch angle ̂ and the output voltage ̂ obtained at the flutter speed of the array for different values of the gap ℎ⁄ are presented in Figure 13. The value of the mean output electric power under the same Based only on the dependence of the flutter speed with the spacing h, two zones, favorable and unfavorable, are identified in Figure 12. It is recalled that each individual harvester in the system only undergoes self-sustaining limit cycle oscillations at or above the flutter speed. Therefore, a favorable characterization is given to the region where lower wind speeds are needed to achieve flutter. A lower flutter critical speed of the system allows the harvesters to operate over a wider Aerospace 2020, 7, 167 20 of 29 range of free-stream wind velocities and to exploit the flutter phenomena at lower wind speeds, something otherwise impossible for the single isolated harvester configuration. Therefore, this scenario is deemed advantageous.
As it is shown in Figure 12, the transition between the favorable and unfavorable interaction zones occurs at a value of spacing h/c approximately equal to five. In the favorable zone, the flutter speed of the array is below the flutter speed of an isolated harvester (V A F < 1). Conversely, in the unfavourable zone, the flutter speed of the array is above the flutter speed of an isolated harvester (V A F < 1). In this last zone, it is observed that a maximum value of V A F is reached approximately at h/c = 7.5. From this value on, on V A F decreases asymptotically approaching the reference value (V A F = 1). This implies that as the distance increases beyond a certain point, the mutual aerodynamic interaction between harvesters becomes insignificant, and each harvester behaves effectively as if it was isolated.
Continuing with the analysis, the amplitudes of the LCO of the wing tip pitch angleθ tip and the output voltagev e obtained at the flutter speed V A F of the array for different values of the gap h/c are presented in Figure 13. The value of the mean output electric power P m under the same conditions is shown Figure 13. These results correspond only to one of the harvesters in the array; thus, the total output power of the system is twice the value presented in Figure 13.
The amplitude of the LCO for the wing tip pitch angle is found to be smaller than the reference value corresponding to an isolated harvester across all the separation range considered. For values of h/c 3, the amplitudeθ tip is observed to decrease with increasing separation until a minimum value at approximately h/c = 3 is reached. From this value on, the amplitude of the pitching motion increases and slowly approaches the reference line.
The values of the amplitude of the output voltagev e (Figure 13b) are well below the baseline value for small values of h/c. For h/c < 2.0, the amplitude of the voltage LCO depicts a steep increase with the separation between harvesters. Approximately for h/c = 2, thev e curve crosses the reference value for an isolated harvester. Beyond this point, the voltage amplitude continues to increase until a maximum value is reached at h/c = 7.5. For h/c > 7.5 the amplitudev e decreases to asymptotically approach the baseline value.
The mean output power (Figure 13c) behaves qualitatively similar to the output voltage LCO amplitude. For values h/c < 3, the presence of the second harvester produces the undesirable outcome of reducing the power harvested by each individual harvester. Approximately for h/c = 3 the mean power P m of each individual harvester in the system is the same as the mean power produced by an isolated harvester. For larger values of h/c, the power continues to grow to reach a maximum at h/c = 7.5 and finally decreases approaching the reference line.
Here, a holistic consideration is elaborated based on the results presented in Figure 12. It is observed that, for separation values h/c < 3, the power output of each individual harvester is below their isolated values. This may suggest that the configuration is detrimental for the energy harvesting performance of each individual harvester. However, such values of the separation h/c lie within what was previously defined as the favorable zone. This means that each harvester is actually undergoing self-sustained LCO (flutter) and generating electrical power at values of the wind speed that would have resulted in no power output at all for a single isolated harvester. Hence, to be fair, a comparison should be made between the zero-power output of a single isolated harvester and the nonzero power output for the array of two harvesters at the given speed. From the above discussion, it is concluded that the definition of the favorable zone answers to the priority given to the ability to harvest energy at lower wind speeds. As a particular example, consider the case of two harvesters separated h/c = 0.5 where the mean power output of each harvester is about 20% lower than the value for an isolated harvester. For the same configuration, the system flutter speed is reduced by 22% with respect to the value for an isolated harvester. Therefore, the array is able to extract energy from a freestream with a 20% smaller speed when the isolated harvester would just simply stay still (or its oscillations damped out). Finally, it is recalled that although the power collected by each individual harvester is Aerospace 2020, 7, 167 21 of 29 20% lower, the combined power extracted by the two-harvester is actually 60% larger than that for a single isolated harvester. This means more power extraction at lower wind speeds but at a lower harvester individual efficiency.
Aerospace 2020, 7, x FOR PEER REVIEW 19 of 28 conditions is shown Figure 13. These results correspond only to one of the harvesters in the array; thus, the total output power of the system is twice the value presented in Figure 13. The amplitude of the LCO for the wing tip pitch angle is found to be smaller than the reference value corresponding to an isolated harvester across all the separation range considered. For values of ℎ⁄ ≲ 3, the amplitude ̂ is observed to decrease with increasing separation until a minimum value at approximately ℎ⁄ = 3 is reached. From this value on, the amplitude of the pitching motion increases and slowly approaches the reference line.

Postcritical Response Analysis
The postcritical behavior of the two-harvester array is analyzed for a fixed value of the dimensionless distance h/c = 5. This separation is chosen because the flutter speed of the system was observed to be approximately equal to the value for an isolated harvester, whereas the voltage and power output are slightly higher due to the mutual aerodynamic interaction. For this particular case, the flutter mode shape, the postflutter behavior, and an analysis of the influence of initial conditions are presented.
As already mentioned, the flutter speed of the two-harvester system with h/c = 5 is equal to the value for one harvester; that is to say, V A F = 41.505 m/s. At this wind speed, because the frequencies of all the modes of the system coalesce into a unique value, the system response can be characterized through what is referred to as the flutter mode shape. This mode shape is presented in Figure 14. It can be observed that the mode consists of coupled bending-torsion motions and that the harvesters are in antiphase (with π radians phase offset), that is to say, the flutter mode is symmetric about the imaginary plane in-between the harvesters. According to the model formulation, the weak aerodynamic interaction between the harvesters contributes sufficient enough nonlinearity to the system to allow postflutter bounded LCO. This is contrary to linear formulations that cannot predict the behavior past the critical speed. For the case under analysis, it is observed that once the flutter speed is exceeded, the system exhibits postcritical limit cycle oscillations. To provide an overall picture of the system response, the bifurcation diagram for the amplitude of the free tip pitch angle is presented in Figure 15. The presence of a supercritical Hopf bifurcation when the wind speed reaches the flutter critical speed is observed in this figure. For wind speeds below the critical speed, the phase diagram of resembles the characteristic response of around a stable focus. In addition, it can be observed that the amplitudes of the LCOs monotonically increase for the range of velocities considered. According to the model formulation, the weak aerodynamic interaction between the harvesters contributes sufficient enough nonlinearity to the system to allow postflutter bounded LCO. This is contrary to linear formulations that cannot predict the behavior past the critical speed. For the case under analysis, it is observed that once the flutter speed is exceeded, the system exhibits postcritical limit cycle oscillations. To provide an overall picture of the system response, the bifurcation diagram for the amplitude of the free tip pitch angle is presented in Figure 15. The presence of a supercritical Hopf bifurcation when the wind speed reaches the flutter critical speed is observed in this figure. For wind speeds below the critical speed, the phase diagram of θ tip resembles the characteristic response of around a stable focus. In addition, it can be observed that the amplitudes of the LCOs monotonically increase for the range of velocities considered.
In Figure 16, a similar bifurcation diagram is presented for the output voltage amplitude. In line with the results for the pitch angle amplitude, the oscillations of the voltage output die away with time when the free-stream speed is below the flutter speed. For wind speeds beyond the flutter speed the output voltage exhibits self-sustained limit cycle oscillations with amplitude monotonically increasing with increasing wind speed.
A final analysis is performed to understand the influence of the initial conditions on the system outputs. For this purpose, the accumulated system electric energy as a function of time t is defined as follows: where P(τ) = v 2 e (τ)/R is the instantaneous combined electrical power of the two harvesters.
limit cycle oscillations. To provide an overall picture of the system response, the bifurcation diagram for the amplitude of the free tip pitch angle is presented in Figure 15. The presence of a supercritical Hopf bifurcation when the wind speed reaches the flutter critical speed is observed in this figure. For wind speeds below the critical speed, the phase diagram of resembles the characteristic response of around a stable focus. In addition, it can be observed that the amplitudes of the LCOs monotonically increase for the range of velocities considered.  In Figure 16, a similar bifurcation diagram is presented for the output voltage amplitude. In line with the results for the pitch angle amplitude, the oscillations of the voltage output die away with time when the free-stream speed is below the flutter speed. For wind speeds beyond the flutter speed the output voltage exhibits self-sustained limit cycle oscillations with amplitude monotonically increasing with increasing wind speed. A final analysis is performed to understand the influence of the initial conditions on the system outputs. For this purpose, the accumulated system electric energy as a function of time is defined as follows: where ( ) = 2 ( )/ is the instantaneous combined electrical power of the two harvesters.
The temporal evolution of is computed for two cases with different set of initial conditions:  Case 1: A 3 deg rotation of the free tip is imposed to both harvesters in a symmetrical configuration. The temporal evolution of E acc is computed for two cases with different set of initial conditions: • Case 1: A 3 deg rotation of the free tip is imposed to both harvesters in a symmetrical configuration. The accumulated energy in time predicted for both cases is shown in Figure 17. It is observed that after a time of approximately 10 s, the accumulated energy in Case 1 is about 47% higher than the corresponding value for Case 2. In addition, the difference between the accumulated energy in both cases grows to reach a constant value of 6.7 J once the LCO amplitude is reached. At this point the curves become parallel to each other.
Aerospace 2020, 7, x FOR PEER REVIEW 23 of 28 as an exchange of energy between the harvesters. The behavior and the longer time to reach the LC may explain why for Case 2 total energy harvested is smaller than for Case 1.
The observed phenomena provide evidence of the relative importance of the initial conditions on the system's harvesting performance when reduced periods of time are considered.  The time histories of the output voltage of each harvester for Cases 1 and 2 are presented in Figure 18. It is observed that for Case 1, the voltage signals of the harvesters are identical (despite the antiphase initial condition) and that both reach the limit cycle (LC) attractor in about 9 s. On the other hand, for Case 2, where only one of the harvesters has initial conditions, it is observed that the LC attractor is reached at a later time compared to the previous case, about 19 s after the start of the motion. Moreover, before reaching the LC, the amplitude of the output voltage of the harvesters appears to be modulated by envelopes that are lagged by 90 deg. This behavior could be interpreted as an exchange of energy between the harvesters. The behavior and the longer time to reach the LC may explain why for Case 2 total energy harvested is smaller than for Case 1.
The observed phenomena provide evidence of the relative importance of the initial conditions on the system's harvesting performance when reduced periods of time are considered. Figure 17. Temporal evolution of the accumulated electric energy of a two-harvester array with separation ℎ⁄ = 5 obtained for two different sets of initial conditions.

Conclusions
This effort is motivated by the need to understand the combined piezoelectric and aeroelastic behavior of arrays of flutter-based energy harvesters. The proposed array configuration is a potential technological solution to overcome the issue of low power outputs of single flutter-based energy harvesting systems. Insights on the effects of the weak aerodynamic coupling between the harvesters as a function of the separation between them are crucial for advancing the field towards reliable, wide operating range and high-power output harvesting systems. The studies presented in this article emphasize the capability of the developed approach to predict the piezoaeroelastic nonlinear response of harvester arrays.
To carry out the present study, a multiphysics framework based on a co-simulation strategy was developed. The three physics involved, fluid dynamics, structural mechanics, and piezoelectricity were represented as two independent, yet interacting, models: a UVLM-based model and computational implementation were developed to compute the aerodynamic forces, and a three-dimensional, a nonhomogeneous piezoelectric/elastic finite-element beam model was introduced to determine the structural response due to the fluid loads. The finite-element beam formulation is based on a classical Euler-Bernoulli beam augmented by the so-called von Kármán nonlinearities. The models were strongly coupled at the computational implementation level through a co-simulation scheme. This approach consists of a bidirectional exchange of information between the physics within an iteration loop at each time step and ensures convergence to a desired precision.
The correct implementation of the computational framework was verified by reproducing a study previously reported by other researchers in the literature. This study consisted of a single flutter-based piezoelectric power-generating wing whose flutter speed and frequency were estimated. The predictions of the current method were shown to agree well with those reported in the literature, thus providing confidence for its use on more complex systems and configurations.
The core study of this effort is on an array of two vertically stacked harvesters. The idea behind this configuration is exploiting the constructive aerodynamic interaction and synergy between the harvesters to enhance the overall performance of the system when compared to the typical single harvester configuration. In this article, the authors focus only on the effect of the vertical separation between harvesters over the cut-in speed of the system (the critical flutter speed at which the system starts harvesting). Other aspects, namely, horizontal separation, number of harvesters, structural coupling between harvesters, end boundary conditions, and so on, will be explored in future works.
From the results obtained, it can be concluded that the cut-in speed of the system can be substantially reduced with respect to the single harvester configuration when the spacing between the vertically stacked harvesters is below a certain threshold. For the particular geometry of the harvesters under study, the separation threshold was found to be approximately five chords. In view of this phenomenon, a favorable interaction region was identified for separation values below this threshold value. An array working in the favorable zone can harvest energy from the wind at speeds that would not induce any postcritical LCO, hence no energy output, on a single harvester alone. In addition, the system can be designed so that the gap between the harvesters can be varied during operation to accommodate the actual wind conditions. Finally, it was also observed that beyond the separation threshold, the cut-in speed and also the power output of each harvester in the array become larger than that of a single harvester configuration. As the separation becomes very large, the aerodynamic coupling between the harvesters becomes very weak, and each device in the array behaves effectively as if they were in a single harvester configuration.
The postcritical behavior of the harvesters is characterized by a supercritical Hopf bifurcation. At postcritical wind speed regimes, the weakly nonlinear aerodynamic interaction between the harvesters was observed to lead to stable LCOs with amplitudes that increase with increasing wind speeds. For these postcritical scenarios, the flutter mode shape could be extracted from the structural response of the system. It was observed that the flutter mode for each of the harvesters is a combination of bending and torsion modes. Moreover, the flutter modes of both harvesters together conform to a symmetric configuration (antiphase).
Finally, the cumulative extracted energy of the system in time was analyzed for two different sets of initial conditions. The goal of this final study was to better understand the transient effects of the nonlinear interaction and the impact of the initial conditions. The results obtained suggest that the initial conditions do have an impact on the elapsed time between the start of the motion and the moment when the stable limit cycle attractor is reached. Based on the observed histories of the accumulated energy, it is concluded that properly set initial conditions may reduce the time needed for the array to reach the attractor for its intended normal operation condition. The closer these initial conditions are to the flutter mode shape, the fastest the LC attractor is reached.
The promising findings reported in this article will help in filling in an important gap regarding the exploitation of constructive aerodynamic interaction and synergies for flutter-based energy harvesting. The results of this effort may find applications in the design and integration of low-energy wind generation technologies into buildings, bridges, and built-in sensor networks in aircraft structures.