## Abstract

Controlling and manipulating elastic/acoustic waves via artificially structured metamaterials, phononic crystals, and metasurfaces have gained an increasing research interest in the last decades. Unlike others, a metasurface is a single layer in the host medium with an array of subwavelength-scaled patterns introducing an abrupt phase shift in the wave propagation path. In this study, an elastic metasurface composed of an array of slender beam resonators is proposed to control the elastic wavefront of low-frequency flexural waves. The phase gradient based on Snell’s law is achieved by tailoring the thickness of thin beam resonators connecting two elastic host media. Through analytical and numerical models, the phase-modulated metasurfaces are designed and verified to accomplish three dynamic wave functions, namely, deflection, non-paraxial propagation, and focusing. An oblique incident wave is also demonstrated to show the versatility of the proposed design for focusing of wave energy incident from multiple directions. Experimentally measured focusing metasurface has nearly three times wave amplification at the designed focal point which validates the design and theoretical models. Furthermore, the focusing metasurface is exploited for low-frequency energy harvesting and the piezoelectric harvester is improved by almost nine times in terms of the harvested power output as compared to the baseline harvester on the pure plate without metasurface.

## 1 Introduction

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*_{0}) 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.

## 2 Theoretical Modeling of the Elastic Metasurface

### 2.1 Introduction to the Metasurface Concept.

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.

*λ*

_{i}and

*λ*

_{t}are the wavelengths of the incident and transmitted waves, respectively, and

*φ*(

*y*) represents the abrupt phase shift introduced by the metasurface. It is necessary to achieve a spatial gradient profile covering the entire 2

*π*phase range to design metasurfaces that can effectively steer the transmitted waves. By properly tailoring each connecting waveguide and satisfying the phase gradient profile, elastic metasurfaces can be designed for the desired wavefronts.

### 2.2 Analytical Formulation of the One-Dimensional Waveguide.

*A*

_{0}) mode Lamb wave propagating in a thin plate with a low Passion ratio is equivalent to the wavenumber of the flexural wave propagating in a beam with the same thickness. Hence, the slender beam resonators (i.e.,

*λ*/

*h*> 10) used in metasurface are modeled with Euler–Bernoulli beam theory, where the effects of the shear deformation and the rotational inertia are neglected. For the unforced vibration, the governing equation of the one-dimensional waveguide can be written as

*w*(

*x*,

*t*) is the vertical displacement in the

*z*-direction,

*EI*is the flexural rigidity,

*ρ*is the density, and

*A*is the area of the cross section. The complete complex solution can be obtained as

*k*= (

*ω*

^{2}

*ρA*/

*EI*)

^{1/4}is the wavenumber of flexural wave,

*ω*is the frequency, and

*j*is the unit imaginary number. The first two terms on the right-hand side of Eq. (3) represent waves propagating in the positive and negative

*x*directions at a phase velocity of

*c*=

*ω*/

*k*. The last two terms represent evanescent waves with exponentially decaying amplitudes with distance. Those waves do not individually transport energy but can be effective and should be considered when the distance between discontinuities (resonator beam length) is significantly smaller than the wavelength as it can impact the metasurface design for maximum wave transmission. The propagating wave (from the left end of the beam) incident on the discontinuity at the boundaries of the thin resonator is scattered, generating various transmitted and reflected waves as depicted in Fig. 2. For an incident wave with unit amplitude, the wave propagation in three regions can be expressed as

*p*and

*e*indicate the propagating and evanescent waves, respectively.

*M*and shear force

*Q*can be written as

*x*= 0 and

*x*=

*L*(where

*L*is the resonator length) and solving the resulting equations simultaneously [44]. Then, the transmission coefficient is calculated by dividing the propagating transmitted wave amplitude by the incident wave amplitude, which contains information about the transmitted phase along with the amplitude as follows:

*T*| is the transmission coefficient and

*ϕ*is the phase of the transmitted wave. The analytical transmission coefficient is shown in Fig. 3(a) for the slender aluminum beam with 50.4 mm length, 13.5 mm width, and 0.794 mm thickness connecting the two aluminum host media with 3.175 mm thicknesses (

*ρ*= 2700 kg/m

^{3},

*E*= 70 GPa). At the frequencies close to the natural frequency of the clamped-clamped beam, the transmission coefficient reaches to unity. In order to verify the transmission coefficient and the analytical modeling framework, we performed a numerical simulation in comsol multiphysics. The time-domain simulations were conducted to obtain the transmission characteristics of the one-dimensional waveguides under sine-burst excitations with different center frequencies. The transmission coefficients were calculated using fast Fourier transformation (FFT) technique and frequency spectrum of the transmission coefficient is obtained as presented in Fig. 3(a). Accordingly, the transmission coefficient results obtained from numerical simulations are in very good agreement with the analytical results. In most engineering applications, including energy harvesting, a higher transmitted coefficient is desired. Considering the highest wave transmission over the metasurface, phase modulation of the slender beams was obtained at 5 kHz by changing the thickness of the beam inserted in the one-dimensional waveguide. As shown in Fig. 3(b), the introduced phase shift covers the entire 2

*π*phase range, which is the fundamental condition for the successful implementation of the proposed metasurface for any desired wave function. We also verified the phase shift of the individual waveguides through numerical simulations as presented with the out-of-plane displacement field in Fig. 3(c). Hence, by properly tailoring beam thicknesses according to the desired phase gradient profile of the metasurface, various anomalous refraction wavefronts can be achieved as explained in the next section.

## 3 Elastic Metasurfaces: Deflecting, Non-Paraxial Propagation, and Focusing Designs

In this section, different metasurface designs are presented along with their verifications to control normally incident *A*_{0} 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.

### 3.1 Wave Deflecting Metasurface.

*A*

_{0}mode Lamb wave, a linear phase gradient in the transverse direction to the wave propagation,

*y*, should be provided. A normal incident propagating wave can be deflected by angle of

*θ*

_{t}when the phase modulation

*φ*(

*y*) is the required phase modulation at location

*y*and is applied as

*φ*(

*y*) is the required phase modulation at location

*y*, and

*k*

_{t}is the wavenumber of the transmitted wave.

Next, we designed the metasurface with 33 slender beams separated at 2 mm gap to deflect normally incident *A*_{0} 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 $\theta 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 $\theta c=\xb1arcsin(2\pi /ktdy\u22121)$, 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*_{0} wave mode at the frequency of 5 kHz.

### 3.2 Non-Paraxial Propagation Metasurface.

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*_{0} 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].

*y*-intercept of the tangent line at any point in the trajectory locate in the metasurface domain. Considering the size of the model we developed, the desired non-paraxial beam trajectory was chosen as the third-order polynomial function fitted to a cubic Bezier curve without loss of generality:

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*_{0} 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.

### 3.3 Wave Focusing Metasurface.

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.

#### 3.3.1 Focusing Metasurface for Waves With Normal Incidence.

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.

#### 3.3.2 Focusing Metasurface for Waves With Oblique Incidence.

*A*

_{0}mode wave at 5 kHz and incident angle $\theta i=20deg$ is considered. In the oblique incident case, the phase modulation property of the waveguide is different than the normal incident case. The detailed derivation is given in the Appendix. In order to focus waves at a focal distance in the transmitted region, the equiphase surface of the transmitted elastic wave should also be a semicircle with the center located at the focus of the metasurface. In the oblique incident case, a correction phase shift should be added to the phase profile which is given as

*φ*

_{c}=

*k*

_{y}

*y*represents the phase retardation of the elastic wave in the incident region, and

*k*

_{y}is the wavenumber component of the incident wave in the

*y*-direction. For the focal distance

*f*set as 152.4 mm, the phase profile and the corresponding thickness profile of the metasurface are presented in Figs. 8(a) and 8(b), respectively.

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.

## 4 Experimental Validation of the Wave Focusing Metasurface

The focusing metasurface was experimentally tested to validate its focusing effect for normally incident *A*_{0} 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.

## 5 Energy Harvesting Metasurface

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*^{2}/*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 $(\u223c1/\omega Cp)$ to harvest normal incident *A*_{0} 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.

## 6 Conclusions

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.

## Acknowledgment

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

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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.

### Appendix: Formulation of Oblique Incidence Flexural Waves

*y*-direction, the general solution can be written as $w(x,y,t)=W(x)ej(\u2212kyy+\omega t)$ where

*k*

_{y}is the wavenumber component of the incident wave in the

*y*-direction. Substituting it into Eq. (A1), it becomes

*k*stands for the wavenumber of the flexural wave in the incident direction.

_{p}*D*=

*Eh*

^{3}/(12(1 −

*ν*

^{2})) represents the bending stiffness of the plate and

*ν*is the Poisson ratio of the material. The definition of other parameters is the same as normally incident case.

*x*-direction. $kx=(kp2\u2212ky2)1/2$ is the wavenumber of the propagating wave in the

*x*-direction. The relationship between

*k*

_{x},

*k*

_{y}, and the angle of incidence,

*θ*

_{i}, is

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

*y*= 0) in different regions can be written as

*p*and

*e*indicate the propagating and evanescent waves, respectively.

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.