Elastic Metasurfaces for Full Wavefront Control and Low-Frequency Energy Harvesting

In the past decade, control and manipulation of acoustic/elastic waves in periodically engineered artificial materials have gained an increasing research attention for their potential in various mechanical, aerospace, and civil engineering systems [1–4]. The bandgap (i.e., ability to forbid wave transmission) [5,6] and dispersive properties (i.e., ability to slow the group velocity) [7] of metamaterials and phononic crystals (PCs) have resulted in unconventional approaches in different applications including sound/vibration attenuation [8,9], cloaking [10–12], filtering [13,14], lensing [15–20], sensing [21,22], subwavelength imaging [23,24], and energy harvesting [25–29]. For instance, Zhu et al. [20] designed a bifunctional PC superlens to simultaneously focus acoustic and flexural waves. Tol et al. exploited gradient-index PC lenses to focus flexural wave energy for enhanced harvesting of the high-frequency flexural wave energy (e.g., tens of kilohertz) [18,19,27]. However, due to the scaling of PCs with the wavelength, such PC-based lens designs would yield very large dimensions to operate at low ambient vibration frequencies. On the other hand, in an effort to harvest low-frequency vibration energy, Chen et al. [28] explored piezoelectric patch array with heavy masses that is based on a PC concept and later Sugino and Erturk [29] investigated a metamaterial harvester composed of piezoelectric cantilevers with tip mass attachments. Alternatively, metasurfaces, which have emerged more recently, offer compact and lightweight solutions in the reduced operating frequencies. Here, our goal is to show the potential of metasurface concepts for full control of waves in plates toward enhanced harvesting of the low-frequency (sub-kHz) elastic wave energy.

A metasurface is composed of an array of subwavelength-scaled structures in the host medium, which can introduce an abrupt phase shift in the wave propagation path and tailor the wavefront based on generalized Snell’s law [30]. By adjusting the phase gradient at the subwavelength scale, elastic/acoustic waves can be manipulated with the smallest possible structure. Since 2014, researchers have studied various acoustic metasurfaces with different structural designs and showed successful manipulation of the acoustic waves [31–34]. On the other hand, elastic metasurfaces are relatively less explored compared to their acoustic counterparts. Zhu and Semperlotti [35] presented the first experimental demonstration of elastic metasurfaces made out of locally resonant torus-like tapers to enable anomalous refraction in solids for both same-mode and mode-converted transmitted Lamb waves at 20.1 kHz. At the same time, Su and Norris [36] proposed an elastic metasurface composed of plate-like waveguides in host medium to control SV- and P-waves in elastic solids. Elastic wave manipulation resulting from bulk wave mode conversion has also been explored by other researchers [37,38]. Moreover, Liu et al. [39] proposed an elastic metasurface composed of zigzag structures to realize source illusion function based on flexural Lamb wave at 8–12 kHz. Lee et al. [40] implemented a sub-structuring design to steer in-plane longitudinal waves at 100 kHz. More recently, Cao et al. [41] presented a pillared elastic metasurface to deflect flexural waves propagating in plates at 6 kHz. Then, they further improved their model to enrich the performance of the metasurface by introducing a disorder concept into their design [42] and by considering an alternative design where there was no slot between the pillared sub-structures [43]. However, most of these studies focused on designing elastic metasurface to control bulk waves and Lamb waves with relatively high frequencies (e.g., tens of kilohertz) or with only a specific wave function (e.g., beam steering).

In this study, we propose an elastic metasurface to fully control and manipulate the elastic wavefront in the sub-kHz regime. The performance of the proposed elastic metasurface is analytically and computationally investigated for wave deflecting, focusing, and non-paraxial propagation of the lowest antisymmetric (A₀) mode Lamb wave. Elastic wave focusing is demonstrated for both normal and oblique incident waves and implemented in enhanced harvesting of low-frequency elastic waves via piezoelectric energy harvesters. The focusing metasurface is fabricated and validated through experimental measurements. Unlike the previous studies that numerically studied the metasurface composed of complex structures, the elastic metasurface proposed in this paper is composed of simple beam resonators. Hence, the analytical solution of the waveguides can be obtained and employed to realize different wave functions with the proposed metasurface concept. This will also enable a deeper understanding of the elastic metasurface mechanism.

The outline of the paper is as follows. In Sec. 2, the analytical model of the one-dimensional waveguide of the metasurface is first developed, and then verified by numerical simulations conducted in comsol multiphysics. Based on this model, the transmission characteristics of the waveguide, including the transmission coefficient and the phase modulation, are obtained for varying thickness of the slender beam. In Sec. 3, based on the phase modulation results, various metasurface designs are proposed to realize different wave functions, including deflecting, focusing, and non-paraxial wave propagation. In addition to the normal incident wave, the focusing design for an oblique incident wave is also presented to demonstrate the versatility of the proposed metasurface. The focusing metasurface is experimentally validated in Sec. 4. Furthermore, in Sec. 5, the focusing design for the normal incident wave is implemented in energy harvesting via piezoelectric energy harvesters to show the potential application of the proposed metasurface. Section 6 summarizes the key points in this study.

The fundamental concept of the metasurface design is to introduce an abrupt phase shift at the interface based on the generalized Snell’s law. As shown in Fig. 1, an array of incident wave impinges on the interface with the angle θ_i. Due to the impedance mismatch at the interface, there are reflected and transmitted waves with angles of θ_r and θ_t, respectively.

In this section, different metasurface designs are presented along with their verifications to control normally incident A₀ mode waves at 5 kHz. To this end, we designed and simulated three metasurfaces realizing different wave functions, namely, wave deflecting, non-paraxial wave propagation, and wave focusing. In addition to normal incident waves, the focusing metasurface concept was also numerically tested for oblique incident waves to show the versatility of the proposed approach.

The proposed metasurface was composed of a finite number of distributed slender beams separated with a constant gap. The thickness of each beam in the metasurface design was properly tailored according to the desired phase gradient and the phase modulation property which was analytically obtained (Fig. 3(b)). Once the metasurface design was finalized, time-dependent simulations were conducted in comsol multiphysics to test the performance of proposed metasurfaces. In the finite element simulations, a high-quality mesh was ensured in the metasurface by setting the mesh size to λ/150. In order to equally resolve the wave in space, a Courant–Friedrichs–Lewy number of 0.2 was selected to set the time-step for the optimal wave solution. In addition, low-reflecting boundary condition was imposed at the edges to reduce the influence of reflected wave component from the plate boundaries.

Next, we designed the metasurface with 33 slender beams separated at 2 mm gap to deflect normally incident A₀ mode wave at 5 kHz by an angle of 20 deg based on the linear phase gradient presented in Fig. 4(a). Figure 4(b) shows the thicknesses of the slender beams in the metasurface that are tailored according to the desired linear phase gradient and the phase modulation property of the waveguide analytically obtained in Sec. 2.

By applying a continuous sinusoidal boundary load at 5 kHz, we numerically tested the metasurface in comsol multiphysics. Figure 4(c) shows the resulting full wavefield, where the transmitted elastic wavefront has been deflected by an angle of $θt=20deg$ as designed. Hence, the time-dependent simulations successfully verify that the metasurface is able to steer the elastic wavefront in the transmitted region by the desired angle.

Note that, due to the phase discretization effect, multiple possible propagating wave solutions can be supported by the same metasurface [35]. The first critical angle for the transmitted wave is $θc=±arcsin(2π/ktdy−1)$⁠, where dy represents the spacing between adjacent beams. Beyond this angle, other refracted beams can be identified in the transmitted region. The number of refracted beams depends on the relationship between dy/λ and the desired angle of the transmitted wave. If dy < λ/2, the critical angle θ_c does not exist. In our design, the spacing between each waveguide, dy, is 16.24 mm and λ is 76.6 mm; hence, there is only one transmitted beam for incident A₀ wave mode at the frequency of 5 kHz.

Another important concept is non-paraxial beam propagation, which can be used as an alternative to the wave cloaking technique. In this section, we present an elastic metasurface for self-bending of the normally incident A₀ wave mode at 5 kHz. To realize non-paraxial propagation in the transmitted region, all of the rays in transmitted region should be tangent to the desired trajectory according to caustic theory [45] as depicted in Fig. 5. For the desired trajectory denoted as a function, y = f(x), the corresponding phase gradient profile φ(y) along the metasurface is obtained with the Legendre transform technique [46].

The required phase gradient profile of the non-paraxial wave metasurface obtained via the Legendre transform technique and its corresponding thickness profile are shown in Figs. 6(a) and 6(b), respectively.

In time-dependent simulations, a continuous sinusoidal force at 5 kHz was vertically applied at the end boundary to excite A₀ mode Lamb wave in the host plate. Figure 6(c) shows the resulting wavefield where the normally incident wave was successfully guided in the non-paraxial beam trajectory in the transmitted region. There is a little discrepancy in the tail of the trajectory due to the discretization effect and the aperture size of the metasurface. The tail trajectory can be further improved with a wider metasurface composed of denser unit cells.

Aside from deflecting and non-paraxial propagation design, the focusing metasurface design could be of particular interest in engineering applications in which localized high intensity wave energy is desired. Hence, we will present a more detailed analysis on the focusing metasurface designs including both normal incidence and oblique incidence waves.

In the focusing metasurface design, we used 35 slender beams constant spacing of 3.5 mm and set the focal distance, f, to 152.4 mm. The required phase gradient profile obtained using Eq. (9) is presented in Fig. 7(a). Then, we designed the metasurface by tailoring the thickness of each waveguide according to their phase modulation properties as shown in Fig. 7(b).

Next, the focusing metasurface was numerically tested by applying a five-cycle sine-burst force with a Gaussian pulse window providing excitations at the center frequency of 5 kHz in finite element simulations. Figure 7(c) shows the instantaneous velocity field where the focusing is clearly seen. Then, the root-mean-square (RMS) velocity around the focal point was calculated by integrating the velocity response over time. Figures 7(d) and 7(e) show the RMS wave field and the RMS velocity normalized by the baseline (i.e., pure plate) simulations along the centerline (y = 0), respectively. Accordingly, the maximum wave intensity occurs at 140.8 mm, with only 7.61% deviation from the theoretical focal point. This error mainly comes from the discretization effect and the aperture size of the metasurface. To reduce this error, a larger metasurface with denser unit waveguide is required. Furthermore, the normalized RMS velocity results show that the vertical wave velocity is magnified by three times in the focal area. Hence, we can verify that the proposed metasurface successfully localizes the sub-kHz range elastic wave energy.

Figures 8(c) and 8(d) show the instantaneous and RMS velocity fields, respectively, verifying the elastic focusing at the desired region. The maximum vertical wave velocity intensity occurs at 154.8 mm, with only a 1.57% derivation from the designed point. Increasing the aperture size of the metasurface designed with denser unit waveguides can further reduce the error. Furthermore, the vertical wave velocity is magnified nearly 3.4 times in the focal area for the $20deg$ angle oblique incident waves as shown in Fig. 8(e). It is very promising that the proposed designs offer significant focusing performance for waves propagating from different incident angles.

The focusing metasurface was experimentally tested to validate its focusing effect for normally incident A₀ mode wave. The metasurface consisting of 35 spatially distributed aluminum beams with width of 13.5 mm, length of 50.4 mm, and spacing of 3.5 mm was fabricated in a 3.175 mm thick aluminum plate by milling process with a machining tolerance of ±0.06 mm. As explained in Sec. 3.3.1, the focal distance of the fabricated metasurface was designed at 152.4 mm. Normally incident plane wave at 5 kHz was generated by exciting 23 equally distributed piezoelectric actuators with length 25 mm, width 5 mm, and thickness 0.3 mm (Steminc) with five cycles of sine-burst signal. A Polytec PSV-500 scanning laser Doppler vibrometer was used to measure the out-of-plane velocity of the focusing region. Note that the overall plate size is chosen large enough to prevent the interference of the transmitted wave packet with the boundary reflections. The experimental measurements were first conducted on the pure plate (without the metasurface) and the baseline results were used in the evaluation of the focusing metasurface. The experimental setup of the plate with focusing metasurface is presented in Fig. 9.

Figures 10(a) and 10(b) show the measured instantaneous and RMS distribution of the out-of-plane velocity fields, demonstrating the focusing phenomenon at the desired region. The maximum vertical wave velocity intensity occurs at 149.0 mm, with only a 2.23% derivation from the designed point. The measured out-of-plane velocity signals at the focus are shown in Fig. 10(c). The experimental normalized RMS velocity results show that the vertical wave velocity is magnified almost 3.1 times in the focal area, demonstrating a good agreement with the numerical simulation results.

One important engineering application with the proposed metasurface is low-frequency wave energy harvesting. To this end, we studied the energy harvester metasurface via multiphysics simulations by modeling a PZT-5A piezoelectric patch (with a thickness of 0.1 mm and half-wavelength diameter of 38.3 mm) and the resistive electrical load (R) and measured the voltage output (v(t)) generated in the piezoelectric harvester as presented in Fig. 11. The electrical power flowing into the resistant was calculated with P = v²/R. We conducted simulations for both metasurface harvester attached at the maximum intensity of the focal region and baseline harvester attached to the pure plate with the optimal resistances $(∼1/ωCp)$ to harvest normal incident A₀ wave excited with sine-burst excitation at 5 kHz.

Figure 11 shows the voltages across the resistor of baseline and metasurface harvesters. It is observed that the harvester voltage is increased by three times using the focusing metasurface design resulting in harvested power increased by 9.04. This significant enhancement in the electrical power is equivalent to the square of normalized RMS velocity as expected.

We proposed an elastic metasurface concept composed of an array of slender beams to control the antisymmetric mode Lamb waves propagating at 5 kHz. The modeling framework for the elastic metasurface was established with both analytical and finite element models. The transmission characteristics and phase modulation property of individual beam resonators in the metasurface were obtained through analytical solution. By properly tailoring the thickness of the beams according to the desired phase gradient profile, three metasurface designs were presented to fully control low-frequency wavefront. To this end, wave deflecting, non-paraxial wave propagation, and wave focusing metasurfaces were numerically tested and verified to accomplish anomalous refraction as desired. Also, the multidirectional wave focusing via the proposed metasurface was demonstrated with an oblique incident case study. Focusing metasurface results showed significant wave focusing with 3 and 3.4 times increase in the wave velocity for the normal and oblique incident cases, respectively. Then, we implemented the proposed metasurface in low-frequency energy harvesting and obtained a dramatic increase in the harvested power output with 9.04 times compared to the baseline harvesting. Aside from energy harvesting, the proposed metasurfaces offer compact solutions in other potential engineering applications at sub-kilohertz frequencies, such as wave circumventing and wireless energy transfer.

This work was supported by the National Science Foundation under Grant No. CMMI-1933436.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtained from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper.

To analyze the phase modulation ability of the metasurface, we can factor out $e−jkyy$ first, and then trace this phase retardation for each waveguide after we obtain the solution to the one-dimensional waveguide located at y = 0.

The unknown complex wave amplitudes in Eq. (A5) are obtained by imposing the linear/angular displacement compatibility and force/moment equilibrium conditions at the two boundaries x = 0 and x = L and simultaneously solving the resulting equations. Then, the transmission coefficient and the phase modulation property can be obtained from $T~3p$ term.