Graded elastic metasurface for enhanced energy harvesting

In elastic wave systems, combining the powerful concepts of resonance and spatial grading within structured surface arrays enable resonant metasurfaces to exhibit broadband wave trapping, mode conversion from surface (Rayleigh) waves to bulk (shear) waves, and spatial frequency selection. Devices built around these concepts allow for precise control of surface waves, often with structures that are subwavelength, and utilise Rainbow trapping that separates the signal spatially by frequency. Rainbow trapping yields large amplifications of displacement at the resonator positions where each frequency component accumulates. We investigate whether this amplification, and the associated control, can be used to create energy harvesting devices; the potential advantages and disadvantages of using graded resonant devices as energy harvesters is considered. We concentrate upon elastic plate models for which the A0 mode dominates, and take advantage of the large displacement amplitudes in graded resonant arrays of rods, to design innovative metasurfaces that trap waves for enhanced piezoelectric energy harvesting. Numerical simulation allows us to identify the advantages of such graded metasurface devices and quantify its efficiency, we also develop accurate models of the phenomena and extend our analysis to that of an elastic half-space and Rayleigh surface waves.


Introduction
Metamaterials, in their modern guise, emerged about two decades ago in optics [1] and their potential in creating artificial media, with properties not found in nature, and the consequent opportunities in wave physics, have triggered intense research activity. Metamaterial concepts have now become a paradigm for the control of waves across much of physics and engineering and are widely used in electromagnetism [1][2][3][4], acoustics [5,6] and elasticity [7]. Whilst many theoretical studies have made the basic concepts well-established, manufacturing complex 3D bulk metamaterials remains a challenge. Attention has therefore focussed on metasurfaces because wave propagation is often dominated by scattering from, transmission through, or waves guided along, surfaces; surfaces are also easier to pattern, print, or manufacture. Metamaterials are available in a variety of shapes and sizes, but all feature arrangements of local heterogeneities are often repeated periodically.
In elasticity, initial designs based around Bragg scattering, phononic crystals, and material contrast were used to create band-gaps [8][9][10] that often drew upon ideas from the photonic crystal community. In geophysical settings, and for applications in Mechanical and Civil engineering, there has been a drive to obtain broader band-gaps at a low frequency [11][12][13][14] and to use locally resonant metamaterials [15][16][17]. Locally resonant metamaterials are characterised by resonating elements, analogous to Helmholtz resonators [18] in acoustics and Fano resonances [19], enabling the creation of low-frequency band-gaps in structures relatively small compared to the wavelength; these ideas can be augmented by tuning or spatially grading the resonators to broaden or decrease band-gaps. At a smaller scale, and consequently higher frequency, wave redirection and protection [20][21][22][23]  Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. tolerances or precision measurements (e.g. interferometry) are required; similar motivations occur in ultrasonic inspection to increase the signal to noise ratio. To capitalize on these recent metamaterial designs, that can improve energy localization over a wide frequency band, energy harvesting is another attractive application; at the millimetric, or centimetric, lengthscales concerned these would operate at frequencies in the kHz range.
Fuelled by demands to reduce the power consumption of small electronic components, vibration-based energy harvesting has received considerable attention over the last two decades; it is attractive to power devices using existing vibrational energy that reduces, or removes, costs associated with periodic battery replacement and the chemical waste of conventional batteries. Many applications arise: wireless sensor networks for civil infrastructure monitoring, unmanned aerial vehicles, battery-free medical sensors implants, and long-term animal tracking sensors [24]. Among the various possible energy harvesting methods, piezoelectric materials offer several advantages due to their large power densities and ease of application [25][26][27]. Piezoelectric energy harvesting is driven by the deformation of the host structure due to mechanical or acoustic vibrations that convert to an electrical potential via embedded piezoelectric materials. To increase harvester efficiency, using ideas based around structuring surfaces, several approaches such as creating a parabolic acoustic mirror, point defects in periodic phononic crystals and acoustic funnels have been employed [28][29][30][31]; lenses to concentrate narrow band vibrations have been proposed using phononic crystals [32][33][34], ideas based around localised defect states [35], and resonant metamaterials endowed with piezoelectric inserts have appeared very recently [36]. Our aim here is to complement these studies by using a graded array to create a metawedge [37], and introduce piezoelectric elements into the array, this addresses one of the main challenges in energy harvesting which is to achieve broadband energy production from ambient vibration spectra. At present, multimodal harvesting, or exploiting nonlinear behaviour, partly addresses this challenge [24]. Multimodal harvesting requires multiple bending modes with close, and effective, resonant peaks thereby leading to a broadband device; these systems can be affected by two problems: low overall power density (power/volume or power/ weight) and complex interface circuits [24]. Similarly to harvesting, resonant metamaterial devices have also been affected by a limited bandwidth. However recent research efforts have proven how this can be enlarged by adopting nonlinear resonators [38,39] and, above all, graded array designs [37].
In this article, our approach is to use graded resonator arrays to concentrate energy exploiting the properties of the Rainbow trapping device shown in figure 1. Because such systems already contain a collection of resonators, the inclusion of vibrational energy harvesters is straightforward leading to truly multifunctional meta-structures combining vibration insulation with harvesting. In section 2 we recall the resonant metawedge design principles for an elastic plate and elastic halfspaces. Section 3 presents the harvester design and we show how it can be tuned by altering the rod length, grading angle and spacing in order to maximise trapping. Here we also introduce a simplified analytical model for the harvester so that optimal design parameters are rapidly found. 3D numerical simulations of the electromechanical device installed on a plate and on a half-space are presented in section 4. Finally, we draw together some concluding remarks in section 5. The elastic metasurface used for energy harvesting (a) and conventional dispersion curves (b) for a periodic array of identical resonators of fixed height (here we take the resonator to have height 581 mm, just before that endowed with piezo-patches in the graded wedge of (a)). Panel (c) illustrates the trapping mechanism for the axial mode due to the metasurface grading by showing the relevant dispersion curve at fixed frequency but varying resonator height; the rod is shown inset at 581mm.

The graded resonant metasurface for elastic waves
A cluster of rods attached to an elastic substrate whether acting as a phononic crystal, for short pillars [40][41][42], or as a resonant metamaterial, for long rods [43], creates a versatile system for controlling elastic waves and lends itself well to the fabrication of graded designs. The physics of the resonant metasurface is described through a Fano-like resonance [19]; a single rod attached to an elastic surface readily couples with the motion of either the A 0 mode in a plate, or the Rayleigh wave on a thick elastic substrate (half-space), and the coupling is particularly strong at the longitudinal resonance frequencies of the rod. Mathematically, the eigenvalues of the equations describing the motion of the substrate and the rod are perturbed by the resonance, mode repulsion occurs, and complex roots arise leading to the formation of a band-gap [44,45]. When the resonators are arranged on a subwavelength cluster (i.e. with λ, the wavelength, much bigger than the resonator spacing), as in the metasurface we use, the resonance of each rod acts constructively enlarging the band-gap until, approximately, the rod's anti-resonance [43,46]. Because the resonance frequency of the rod, which depends on rod height, determines the band-gap position, then by simply varying the length of the rods one gets an effective band-gap that is both broad and subwavelength; the rod length is therefore a key parameter for metasurface tunability. However both periodicity and/or the distribution of the rods cannot be overlooked as they also shape the dispersion curves leading to a zone characterised by dynamic anisotropy and negative refraction [47,48]. The band-gap frequency f is calculated approximately starting from the resonator length, h, and the well-known formula [20,48]: where E is the Young's modulus of the material and ρ its density. Equation (1) holds provided the substrate is sufficiently firm so that one end of the rod is effectively clamped. Conversely, an overly soft substrate (e.g. a very thin plate) may result in a rod behaving as a rigid body connected with a spring thus compromising the amplification. In practical terms, this provides bounds on the plate thickness and requires cross-checking, via finite element simulations, the modal shape of the longitudinal eigenmode. Similar considerations are valid for the half-space, that, again, cannot be too soft when compared to the resonators [46]. A schematic showing the resonant metawedge is depicted in figure 1(a). We place a graded array of rods atop an elastic plate, and these rods have resonances (both flexural and longitudinal); if we assume the rods have constant height then these naturally appear in the dispersion curves of figure 1(b). These dispersion curves are computed along the 2D irreducible Brillouin zone [49] using the finite elements code COMSOL Multiphysics®, in a physical cell containing a single resonator, that incorporate the Bloch phase shift via Bloch-Floquet periodic boundary conditions. Figure 1(b) shows the longitudinal resonance (dashed line), and additionally a large number of flexural resonances are also visible. The impact of these flexural resonances on the wave propagation is negligible [50], at least in terms of the band-gaps. An unconventional dispersion curve representation, where the frequency is fixed and the resonator height varied, as in figure 1(c), shows how rod height affects the array. By grading the array from short to tall resonators, figure 1(c), we see that the dispersion curve hits the bottom of the band-gap and subsequently the energy must be reversed and propagate back through the array; equally importantly the group velocity tends to zero at this turning point and the wave slows down and therefore spends time in the vicinity of the resonator feeding energy into it.
The resonant metawedge can also be attached onto a deep elastic substrate where it is capable of mode conversion or the trapping and reflection of Rayleigh waves. The effect is underpinned by the hybridisation between the longitudinal resonance of the rods, and the vertical component of Rayleigh waves [37,43,46]; the type of interaction, conversion or trapping, depends on the direction of wave propagation across the metawedge. When moving towards increasing rod lengths we obtain trapping, conversely by reversing the wave direction, we obtain Rayleigh to shear, S, wave conversion. Because each rod has a different resonant frequency, the turning point (where conversion or trapping take place) varies spatially according to the frequency content of the incident wave (the so-called Rainbow effect [51,52]). Trapping is by far the most interesting property for energy harvesting as it leads to higher amplification on the resonator responsible for the trapping. For the elastic plate, as the substrate, there is no conversion to a bulk mode as this requires a deep substrate, but the concept of trapping carries across to the antisymmetric A 0 Lamb mode.

Harvester design
Concentrating energy at a known spatial position, as the graded array does, is only part of the requirement for harvesting: it is also necessary to design an arrangement for the piezoelectric patches that takes full advantage of the displacements induced within the array and upon the rods. To understand the mechanisms involved we will take a single piezo-augmented resonator and place it within the graded array (figure 1).
Recalling that the dominant mode of interest is the longitudinal one, the harvester we use is the rod, with four cantilever beams arranged in a cross-like shape placed upon the top, as depicted in figure 2(a); each beam embeds a piezoelectric patch and their motion in harvesting mode is shown in figure 2(b). In our elastic modelling we consider the plate and rods to be made from aluminium (E a =70 GPa, ν a =0.33 and ρ a =2710 kg m −3 ) while the piezoelectric patches are made of PZT-5H (E p =61 GPa, ν p =0.31 and ρ p =7800 kg m −3 ). The harvester design, figure 2(a), has strong dynamic coupling between the rod and the cantilever beams as we carefully design it to work in the double amplification regime coupling the axial fundamental mode with the flexural one of the beam (see figure 2(b)) i.e. the lengths of the cantilever beams are not arbitrary. Figure 2(c) illustrates the transmission spectra obtained from a simple spring-mass model (see inset). Three transmission spectra are shown depending on the output and input position: harvester (w b /w 0 ), cantilever beam (w b /w r ) and rod with beams (w r /w 0 ). The resonance of the cantilever beam coincides with the antiresonance of the rod with beams. At this frequency, because of energy conservation, the energy is localised in the cantilever beam.
Provided the axial and flexural frequencies are matched, an amplification over a wide frequency range is obtained. We design the harvester to maximize this interaction at the trapping frequency of the metasurface in a range between 2 and 2.4 kHz. We use an optimization procedure that maximizes the integral of the transmission spectrum in this frequency range subject to an area constraint given by the geometry of the piezoelectric patch (defined on the basis of commercial values) whose width and length are the same as the beam. The values are: A r =25 mm 2 (square cross section), A b =10 mm 2 (same width as the rod, i.e. 5 mm), l r =460 mm, l b =25 mm and h p =0.3 mm (the subscript r denotes the properties of the rod, b of the beams and p of the patches). This harvester design provides a multimodal response of the system, always exploiting the fundamental flexural mode of the beams. The fundamental mode of the cantilever beam is that most suited for energy harvesting as it does not result in charge cancellation and has the highest power production; with higher mode numbers this would decay by around two orders of magnitude [24,25]. In addition, this configuration provides high amplification at low frequency (relevant for ambient vibrations) due to both the slenderness of the rod and the added mass from the four cantilever beams. Tuning the cantilever beams to couple optimally with the fundamental mode of the rod leads to a rather lengthy numerical exercise, this is side-stepped here by using a simple, yet highly effective, mass-spring model. We outline the model below and note its accurate predictions in figure 3 in the frequency range of interest.
We use the 1D spring-mass model, inset in figure 2(c), and adopt Timoshenko beam theory. Defining with w r (z, t) and w b (x, t) the vertical displacement of the rod and the beam respectively, the effective mass is obtained writing the kinetic energy of the system for both the axial and flexural mode as: while elastic and electric lumped coefficients directly follow from the internal energy definition: with a separation of variables taken as: ( ) ( ) = w z t w t z l , p p the total mass of the four beams, T and S the stress and strain fields, θ the electromechanical coupling coefficient and C p the internal piezoelectric capacitance. Lumped coefficients are then defined as: a a a being the aluminium shear modulus, a p the patch width, b b the beam thickness, y n the neutral axis of the composite beam, e 31 the considered piezoelectric coefficient and 33 the constant-stress dielectric constant. The dynamics, of the electromechanical problem, are now succinctly described by three linear coupled ordinary differential equations: ⎧ where R is the electrical resistance. Fourier transforming (6) gives the transfer function as: from which the voltage and power follow directly as: with w 0 the imposed harmonic displacement amplitude. For brevity we do not incorporate damping, but it can be easily included in (6) by introducing a term linearly dependent on the velocity  w. The main contributor to damping in this system is provided by the piezoelectric material since the quality factor of aluminium is very high. However we neglected damping here since the thickness of the patch is very small with respect to that of the beam. This simplified model provides surprisingly accurate predictions in the frequency range of interest (see figures 3(a) and (b)), with just a 0.7% of error in the natural frequency prediction. The product RC p defines the time constant of the circuit τ RC providing measure of the time required to charge the capacitor through the resistor. The time constant is related to the cutoff circular frequency which is an alternative parameter of the RC circuit: The values of R maximising the electric power at each excitation frequency are obtained by imposing its stationarity with respect to R (the dashed white line in figure 4(a)): opt. The optimal electrical resistance is then obtained from the intersection of (8) and (9) (the intersection of black and white dashed lines in figure 4 (a)) i.e. = W R 7.08k opt . Now that we have obtained a harvester optimised for longitudinal rod resonance, it can be introduced into the metasurface array to assess whether its performance is increased when Rainbow trapping occurs. The considerable prior research with this metamaterial in the frequency range of interest [21,43] provides essential insights to the numerical simulations.

Numerical results
We present numerical simulations of the graded harvester on an elastic plate, the plate has the form of a finite width strip to reduce computational complexity [37,50] whilst still preserving the physics of the resonant metasurface. The array is attached to an aluminium strip 30mm wide and of thickness of 10 mm which is sufficiently stiff to avoid anomalous resonances in the rod. The array is composed of 30 rods, each with square cross section of area 25 mm 2 , a linear height gradient ranging from 250 to 650 mm and a constant spacing between the rods of 15 mm; this results in a ∼43°slope angle. Equation (1) suggests that the metasurface will have the longitudinal fundamental mode in the range 2 to 5 kHz and rod number 26 (with the rod numbering in the array having 1 the shortest resonator and 30 the longest), with length 460 mm, is the harvester designed in section 3.
To benchmark the graded harvester, and to highlight the advantage of the elastic wave trapping, a time domain simulation with a Gaussian excitation centered at 2.15 kHz (frequency of the main mode of the harvester) is performed. This is compared to an isolated harvester on the same plate with, and without, the metasurface (see figures 5(a) and (b)) and to the conventional harvesting scheme where the cantilever beams are directly mounted on the waveguide (see inset (c) in figure 5(c)).
Both models are excited with a vertical pressure p 0 , on a surface S 0 =650 mm 2 , applied in the same position which generates an A 0 Lamb wave. The amplitude of the pressure is such that a maximum acceleration of 1g is imposed at the input. FEM simulations on these 3D structures are performed in ABAQUS CAE 2018 introducing piezoelectric elements on the top of the harvester, with zero voltage boundary conditions on the bottom face of the piezoelectric patches. Electrodes are modelled introducing a voltage constraint on the piezofaces in order to define equipotential surfaces. The electrodes on top are connected in series and attached to the previously defined optimal resistance; to model this, a customised Fortran subroutine has been implemented in ABAQUS, according to [53,54]. Spurious edge reflections are avoided by imposing absorbing boundary conditions at the end of the strip along x direction. They are implemented in ABAQUS using the ALID (Absorbing Layers using Increasing Damping) method, adopting a cubic law function for mass proportional Rayleigh damping with = c 10 M max 6 , as described in [55] and as widely used in the elastic wave community. It can be clearly seen that the metasurface (see figure 5(c)) provides a strong amplification of the electrical power, and when the harvester is embedded in the metasurface, due to the lower group velocity of the waves interacting with the array grading, the power generation peak occurs at a later time. The wave velocity decreases along the grading, reaching the harvester with a nearly zero value. For short propagation time (e.g. 10 ms), the grading has a detrimental effect on the power production. Conversely if the excitation is sufficiently long, energy stops in the position of the harvester loading it continuously and, in this case, the energy production peaks at about 20 ms; after this time, waves are backscattered and the power production decreases. It is worth noticing that not only is the peak higher, but globally the area below the curve is larger: this further demonstrates that energy is trapped inside the metasurface; trapping strongly depends on the grading angle [56]. Having a very small angle means a lower variation of the axial resonance frequency between each rod. On one hand this leads to a smooth group velocity variation reducing possible scattering, on the other it has a negative effect on the harvestable energy bandwidth as the rod resonant peaks will likely overlap. Another trade off must be found between having a sufficient number of rods before the harvester, to smoothly reduce the wave velocity while, at the same time, minimising the leakage due to parasite flexural rods vibrations (e.g. figure 1(b)) that reduce the energy available to the harvester(s). We recall that the metasurface also provides strong vibration attenuation over the entire frequency range due to bandgap generated by the longitudinal resonances of the rods. The behaviour is strongly subwavelength, as can be noticed comparing the rod spacing with the wavelength of the A 0 mode in the frequency range we consider.
It is possible to increase the harvesting bandwidth by introducing other harvesters in different positions. In this way, for broadband input, energy is almost uniformly trapped inside the grading and feeds different harvesters. The harvesters introduced maintain the same beams on top and only differ in terms of the rod length. From a practical, and economical perspective, this solution is more feasible since the piezo patches are equal and so is the optimal resistance R opt . Theoretically, a beam grading induces higher power production for the isolated harvesters, as explained in section 3. However, since the harvesters are electrically connected in series, the probability of out of phase responses and hence charge cancellation increases. Therefore, unless one aims at powering multiple devices accommodating different R opt , keeping the same piezoelectric beam size leads to a simpler and more controllable system. The rod length of the harvesters is defined adopting a linear variation starting from resonator 26, with the same metasurface angle. As for the previous case, the piezoelectric patches of each harvester are connected in series to the same electrical resistance. The input vibration is a zero mean coloured noise with a frequency content mainly in the metasurface range between 2 and 5 kHz. This is applied with a uniform pressure on the same surface S 0 , such that a maximum 1g acceleration is imposed at the input. The optimal electrical resistance in this case is obtained with a parametric analysis of the entire system. This means that (9), derived for an harmonic regime, cannot be used here because of the dynamic interaction between harvester of different length, the broadband and transient nature of the input signal; by performing a grid search over several simulations an optimal value of 1.18 kΩ was obtained. The case with isolated harvesters is then compared with the one where they are introduced inside the metasurface (figures 6(a) and (b)). As in the previous case, it can be clearly seen that the metasurface (figure 6(c)) strongly amplifies the electrical power. The global behaviour shown in figure 6(c) is really similar to the one reported in figure 5(c). The longer time necessary to reach the power production peak points to the same energy trapping mechanism discussed for the single harvester. Finally, we remark a phase synchronisation across the harvesters enhanced by the presence of the grading that partially avoids charge cancellation when they are connected in series.
The same mechanism is found when placing the same metasurface on an elastic half-space. As for the plate we took a finite width strip, but now of infinite depth with plane strain boundary conditions applied on these faces along the y direction. This set-up represents the case where the metasurface is built directly on the ground or on a thick mechanical component to harvest ambient noise. We use the same coloured noise excitation adopted for the plate strip, adopting a uniform pressure able to impose a maximum 1g acceleration at the input. The same harvesters are placed inside, and outside, the metasurface (see figures 7(a) and (b)) and connected in series with an electrical resistance equal to the previous one (it is assumed to power the same device). The metasurface is able to provide a strong localization due to Rainbow trapping, inducing, at the same time, a Rayleigh wave bandgap (figures 7(a) and (b)). The different level of amplification provided can be attributed to the much lower unit-cell size to wavelength ratio for Rayleigh waves as well as the strong energy scattering in the substrate that for the plate did not occur.

Concluding remarks
We have numerically demonstrated the advantage of combining graded resonant metamaterials with multiple piezoelectric harvesters both in terms of band widening and overall efficiency. This is demonstrated for both a plate, and a half-space, with a graded array resonators; this exploits the wave trapping mechanism, increasing the performance of the introduced piezoelectric energy harvesters. To exploit this behaviour we find that a sufficiently long excitation is required as the wave group velocity decreases along the grading reaching the value of zero in the position of the harvester. In this manner, before being backscattered, the waves enjoy a longer interaction with the harvester enhancing its power production when compared against the case of isolated harvesters. Additional features are that the system can be frequency-tuned simply by adding masses on the top of the rods and in this way it is possible to match the longitudinal mode of the rod and the flexural one of the piezoelectric patch even at very low frequencies; this is important from an application point of view where one wish to exploit actual low frequency ambient spectra (much lower than 1 kHz). Future work will act to address both the experimental verification of such a device, the optimisation of the graded profile and the integration of nonlinear elements in the array or on the plate capable to exploit even lower frequency vibrations.