## 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 [14]. 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 [1012], filtering [13,14], lensing [1520], sensing [21,22], subwavelength imaging [23,24], and energy harvesting [2529]. 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 [3134]. 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 (A0) 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.

Fig. 1
Fig. 1
Close modal
The generalized Snell’s law can be written as [30]
$12πdφ(y)dy=1λtsinθt−1λisinθi$
(1)
where λ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.

The metasurface proposed in this study is composed of an array of slender beam resonators, which act like phase modulators and are used to connect two infinite plates, as illustrated in Fig. 2. At low frequencies, the wavenumber of the lowest antisymmetric (A0) 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
$∂4w(x,t)∂x4+ρAEI∂2w(x,t)∂t2=0$
(2)
where 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
$w(x,t)=(A~e−jkx+B~ejkx+C~e−kx+D~ekx)ejωt$
(3)
where 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
$w(x,t)={(I~e−jk1x+R~1pejk1x+R~1eek1x)ejωt,x∈(−∞,0)(T~2pe−jk2x+T~2ee−k2x+R~2pejk2x+R~2eek2x)ejωt,x∈(0,L)(T~3pe−jk3x+T~3ee−k3x)ejωt,x∈(L,∞)$
(4)
where $R~$ and $T~$ are the complex amplitudes of the reflected and transmitted wave components, respectively, and $I~$ is the complex amplitude of the incident wave. The subscripts 1, 2, and 3 represent the wave components propagating in the corresponding regions, and the subscripts p and e indicate the propagating and evanescent waves, respectively.
Fig. 2
Fig. 2
Close modal
The bending moment M and shear force Q can be written as
$M=−EI∂2w∂x2,Q=−EI∂3w∂x3$
(5)
The complex wave amplitudes are obtained by imposing the linear/angular displacement compatibility and force/moment equilibrium conditions at the resonator boundaries at 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=|T|ejϕ$
(6)
where |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/m3, 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.
Fig. 3
Fig. 3
Close modal

## 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 A0 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.

According to the generalized Snell’s law, in order to deflect the incoming A0 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)=y×ktsinθt$
(7)
where φ(y) is the required phase modulation at location y, and kt 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 A0 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.

Fig. 4
Fig. 4
Close modal

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 A0 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 A0 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].

Fig. 5
Fig. 5
Close modal
In order to obtain a high-resolution non-paraxial beam in the transmitted region, the proposed metasurface should be long enough to ensure the 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:
$y=f(x)=1.59x3−3.43x2+1.67x−0.25$
(8)

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.

Fig. 6
Fig. 6
Close modal

In time-dependent simulations, a continuous sinusoidal force at 5 kHz was vertically applied at the end boundary to excite A0 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.

To achieve wave focusing at a focal distance in the transmitted region, the equiphase surface of the transmitted wave should be a semicircle with the center located at the focus. Hence, all rays emitted from the metasurface have the same phase when they arrive at the focus. To obtain such an equiphase surface in the transmitted region, we need to tailor the phase profile into a hyperbolic form [35] given as
$φ(y)=kt(y2+f2−f)$
(9)

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).

Fig. 7
Fig. 7
Close modal

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.

The lowest antisymmetric mode Lamb wave with oblique incidence is also studied to demonstrate the versatility of the proposed metasurface in controlling incident waves propagating from different incident angles. In this section, the oblique incident A0 mode wave at 5 kHz and incident angle $θ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
$φ(y)=kt(y2+f2−f)+φc$
(10)
where φc = kyy represents the phase retardation of the elastic wave in the incident region, and ky 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.
Fig. 8
Fig. 8
Close modal

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 A0 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.

Fig. 9
Fig. 9
Close modal

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.

Fig. 10
Fig. 10
Close modal

## 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 = v2/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 A0 wave excited with sine-burst excitation at 5 kHz.

Fig. 11
Fig. 11
Close modal

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

The analytical model of the waveguide of the oblique incident case can be established using Kirchhoff plate theory [47]. The governing equation of the Kirchhoff plate can be written as
$∂4w∂x4+2∂4w∂x2∂y2−∂4w∂y4+ρhD∂2w∂t2=0$
(A1)
For oblique incident waves, due to the phase-matching property in the y-direction, the general solution can be written as $w(x,y,t)=W(x)ej(−kyy+ωt)$ where ky is the wavenumber component of the incident wave in the y-direction. Substituting it into Eq. (A1), it becomes
$d4Wdx4−2ky2d2Wdx2−η4W=0$
(A2)
where $η4=kp4−ky4$ and $kp4=ω2ρh/D$. Here, kp stands for the wavenumber of the flexural wave in the incident direction. D = Eh3/(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.
The solution to Eq. (A1) can be written as
$w(x,y,t)=(A~ej(−kyy−kxx)+B~ej(−kyy+kxx)+C~e−kxexe−jkyy+D~ekxexe−jkyy)ejωt$
(A3)
where $kxe=(kp2+ky2)1/2$ stands for the wavenumber of the evanecence wave in the x-direction. $kx=(kp2−ky2)1/2$ is the wavenumber of the propagating wave in the x-direction. The relationship between kx, ky, and the angle of incidence, θi, is
$kxcosθi=kysinθi$
(A4)

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.

Assuming a unit amplitude of the vertical displacement of the incident wave, the wave solutions for the middle waveguide (y = 0) in different regions can be written as
$w(x,y,t)={(e−jkxx+R~1pejkxx+R~1eekxex)ejωt,x∈(−∞,0)(T~2pe−jk2x+T~2ee−k2x+R~2pejk2x+R~2eek2x)ejωt,x∈(0,L)(T~3pe−jkxx+T~3ee−kxex)ejωt,x∈(L,∞)$
(A5)
where $T~$ represents transmitted wave, and $R~$ represents reflected wave. The subscripts p and e indicate the propagating and evanescent waves, respectively.
The moments and shear force of the plate per unit width can be defined as
$Mx=∫−h/2h/2σxzdz=−D(∂2w∂x2+ν∂2w∂y2)$
(A6)
$Mxy=−∫−h/2h/2τxyzdz=D(1−ν)∂2w∂x∂y$
(A7)
$Qx=−D∂∂x(∂2w∂x2+∂2w∂y2)$
(A8)
The net vertical force per unit width can be written as
$Vx=Qx−∂Mxy∂y$
(A9)

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.

## References

1.
Lu
,
M. H.
,
Feng
,
L.
, and
Chen
,
Y. F.
,
2009
, “
Phononic Crystals and Acoustic Metamaterials
,”
Mater. Today
,
12
(
12
), pp.
34
42
. 10.1016/S1369-7021(09)70315-3
2.
Hussein
,
M. I.
,
Leamy
,
M. J.
, and
Ruzzene
,
M.
,
2014
, “
Dynamics of Phononic Materials and Structures: Historical Origins, Recent Progress, and Future Outlook
,”
ASME Appl. Mech. Rev.
,
66
(
4
), p.
040802
. 10.1115/1.4026911
3.
Ma
,
G.
, and
Sheng
,
P.
,
2016
, “
Acoustic Metamaterials: From Local Resonances to Broad Horizons
,”
,
2
(
2
), p.
e1501595
4.
Jin
,
Y.
,
Djafari-Rouhani
,
B.
, and
Torrent
,
D.
,
2019
, “
Gradient Index Phononic Crystals and Metamaterials
,”
Nanophotonics
,
8
(
5
), pp.
685
701
. 10.1515/nanoph-2018-0227
5.
Deymier
,
P. A.
,
2013
,
Acoustic Metamaterials and Phononic Crystals
,
Springer-Verlag
,
Berlin
.
6.
Laude
,
V.
,
2015
,
Phononic Crystals: Artificial Crystals for Sonic, Acoustic, and Elastic Waves
,
Walter de Gruyter GmbH and Co KG.
,
Berlin
.
7.
Page
,
J. H.
,
Sheng
,
P.
,
Schriemer
,
H. P.
,
Jones
,
I.
,
Jing
,
X.
, and
Weitz
,
D. A.
,
1996
, “
Group Velocity in Strongly Scattering Media
,”
Science
,
271
(
5249
), pp.
634
637
. 10.1126/science.271.5249.634
8.
Matlack
,
K. H.
,
Bauhofer
,
A.
,
Krödel
,
S.
,
Palermo
,
A.
, and
Daraio
,
C.
,
2016
, “
Composite 3D-Printed Metastructures for Low-Frequency and Broadband Vibration Absorption
,”
,
113
(
30
), pp.
8386
8390
. 10.1073/pnas.1600171113
9.
Li
,
Y.
,
Baker
,
E.
,
Reissman
,
T.
,
Sun
,
C.
, and
Liu
,
W. K.
,
2017
, “
Design of Mechanical Metamaterials for Simultaneous Vibration Isolation and Energy Harvesting
,”
Appl. Phys. Lett.
,
111
(
25
), p.
251903
. 10.1063/1.5008674
10.
Popa
,
B. I.
,
Zigoneanu
,
L.
, and
Cummer
,
S. A.
,
2011
, “
Experimental Acoustic Ground Cloak in Air
,”
Phys. Rev. Lett.
,
106
(
25
), p.
253901
. 10.1103/PhysRevLett.106.253901
11.
Craster
,
R. V.
, and
Guenneau
,
S.
,
2012
,
Acoustic Metamaterials: Negative Refraction, Imaging, Lensing and Cloaking
,
Springer
,
London
.
12.
Lee
,
M. K.
, and
Kim
,
Y. Y.
,
2016
, “
Add-On Unidirectional Elastic Metamaterial Plate Cloak
,”
Sci. Rep.
,
6
(
1
), pp.
1
10
.
13.
Pennec
,
Y.
,
Djafari-Rouhani
,
B.
,
Vasseur
,
J. O.
,
Khelif
,
A.
, and
Deymier
,
P. A.
,
2004
, “
Tunable Filtering and Demultiplexing in Phononic Crystals With Hollow Cylinders
,”
Phys. Rev. E
,
69
(
4
), p.
046608
. 10.1103/PhysRevE.69.046608
14.
Lee
,
H. J.
,
Lee
,
J. K.
, and
Kim
,
Y. Y.
,
2015
, “
Elastic Metamaterial-Based Impedance-Varying Phononic Bandgap Structures for Bandpass Filters
,”
J. Sound Vib.
,
353
, pp.
58
74
. 10.1016/j.jsv.2015.05.012
15.
Lin
,
S. C. S.
,
Huang
,
T. J.
,
Sun
,
J. H.
, and
Wu
,
T. T.
,
2009
, “
,”
Phys. Rev. B
,
79
(
9
), p.
094302
. 10.1103/PhysRevB.79.094302
16.
Climente
,
A.
,
Torrent
,
D.
, and
Sánchez-Dehesa
,
J.
,
2010
, “
Sound Focusing by Gradient Index Sonic Lenses
,”
Appl. Phys. Lett.
,
97
(
10
), p.
104103
. 10.1063/1.3488349
17.
Jin
,
Y.
,
Torrent
,
D.
,
Pennec
,
Y.
,
Pan
,
Y.
, and
Djafari-Rouhani
,
B.
,
2016
, “
Gradient Index Devices for the Full Control of Elastic Waves in Plates
,”
Sci. Rep.
,
6
(
1
), p.
24437
. 10.1038/srep24437
18.
Tol
,
S.
,
Degertekin
,
F. L.
, and
Erturk
,
A.
,
2017
, “
Phononic Crystal Luneburg Lens for Omnidirectional Elastic Wave Focusing and Energy Harvesting
,”
Appl. Phys. Lett.
,
111
(
1
), p.
013503
. 10.1063/1.4991684
19.
Tol
,
S.
,
Degertekin
,
F. L.
, and
Erturk
,
A.
,
2019
, “
3D-Printed Phononic Crystal Lens for Elastic Wave Focusing and Energy Harvesting
,”
,
29
, p.
100780
.
20.
Zhu
,
Y.
,
Cao
,
L.
,
Merkel
,
A.
,
Fan
,
S. W.
, and
Assouar
,
B.
,
2020
, “
Bifunctional Superlens for Simultaneous Flexural and Acoustic Wave Superfocusing
,”
Appl. Phys. Lett.
,
116
(
25
), p.
253502
. 10.1063/5.0004428
21.
Chen
,
Z.
,
Guo
,
B.
,
Yang
,
Y.
, and
Cheng
,
C.
,
2014
, “
Metamaterials-Based Enhanced Energy Harvesting: A Review
,”
Phys. B
,
438
, pp.
1
8
. 10.1016/j.physb.2013.12.040
22.
Danawe
,
H.
,
Okudan
,
G.
,
Ozevin
,
D.
, and
Tol
,
S.
,
2020
, “
Conformal Gradient-Index Phononic Crystal Lens for Ultrasonic Wave Focusing in Pipe-Like Structures
,”
Appl. Phys. Lett.
,
117
(
2
), p.
021906
. 10.1063/5.0012316
23.
Sukhovich
,
A.
,
Merheb
,
B.
,
Muralidharan
,
K.
,
Vasseur
,
J. O.
,
Pennec
,
Y.
,
Deymier
,
P. A.
, and
Page
,
J. H.
,
2009
, “
Experimental and Theoretical Evidence for Subwavelength Imaging in Phononic Crystals
,”
Phys. Rev. Lett.
,
102
(
15
), p.
154301
. 10.1103/PhysRevLett.102.154301
24.
Zhu
,
J.
,
Christensen
,
J.
,
Jung
,
J.
,
Martin-Moreno
,
L.
,
Yin
,
X.
,
Fok
,
L.
,
Zhang
,
X.
, and
Garcia-Vidal
,
F. J.
,
2010
, “
A Holey-Structured Metamaterial for Acoustic Deep-Subwavelength Imaging
,”
Nat. Phys.
,
7
(
1
), pp.
52
55
. 10.1038/nphys1804
25.
Gonella
,
S.
,
To
,
A. C.
, and
Liu
,
W. K.
,
2009
, “
Interplay Between Phononic Bandgaps and Piezoelectric Microstructures for Energy Harvesting
,”
J. Mech. Phys. Solids
,
57
(
3
), pp.
621
633
. 10.1016/j.jmps.2008.11.002
26.
Chen
,
T.
,
Li
,
S.
, and
Sun
,
H.
,
2012
, “
Metamaterials Application in Sensing
,”
Sensors
,
12
(
3
), pp.
2742
2765
. 10.3390/s120302742
27.
Tol
,
S.
,
Degertekin
,
F. L.
, and
Erturk
,
A.
,
2016
, “
Gradient-Index Phononic Crystal Lens-Based Enhancement of Elastic Wave Energy Harvesting
,”
Appl. Phys. Lett.
,
109
(
6
), p.
063902
. 10.1063/1.4960792
28.
Chen
,
Z.
,
Yang
,
Y.
,
Lu
,
Z.
, and
Luo
,
Y.
,
2013
, “
Broadband Characteristics of Vibration Energy Harvesting Using One-Dimensional Phononic Piezoelectric Cantilever Beams
,”
Phys. B
,
410
, pp.
5
12
. 10.1016/j.physb.2012.10.029
29.
Sugino
,
C.
, and
Erturk
,
A.
,
2018
, “
Analysis of Multifunctional Piezoelectric Metastructures for Low-Frequency Bandgap Formation and Energy Harvesting
,”
J. Phys. D: Appl. Phys.
,
51
(
21
), p.
215103
. 10.1088/1361-6463/aab97e
30.
Yu
,
N.
,
Genevet
,
P.
,
Kats
,
M. A.
,
Aieta
,
F.
,
Tetienne
,
J.-P.
,
Capasso
,
F.
, and
Gaburro
,
Z.
,
2011
, “
Light Propagation With Phase Discontinuities: Generalized Laws of Reflection and Refraction
,”
Science
,
334
(
6054
), pp.
333
337
. 10.1126/science.1210713
31.
Assouar
,
B.
,
Liang
,
B.
,
Wu
,
Y.
,
Li
,
Y.
,
Cheng
,
J. C.
, and
Jing
,
Y.
,
2018
, “
Acoustic Metasurfaces
,”
Nat. Rev. Mater.
,
3
(
12
), pp.
460
472
. 10.1038/s41578-018-0061-4
32.
Xie
,
Y.
,
Wang
,
W.
,
Chen
,
H.
,
Konneker
,
A.
,
Popa
,
B. I.
, and
Cummer
,
S. A.
,
2014
, “
Wavefront Modulation and Subwavelength Diffractive Acoustics With an Acoustic Metasurface
,”
Nat. Commun.
,
5
(
1
), pp.
1
5
.
33.
Li
,
Y.
,
Jiang
,
X.
,
Liang
,
B.
,
Cheng
,
J.-C.
, and
Zhang
,
L.
,
2015
, “
Metascreen-Based Acoustic Passive Phased Array
,”
Phys. Rev. Appl.
,
4
(
2
), p.
024003
. 10.1103/PhysRevApplied.4.024003
34.
Ma
,
G.
,
Yang
,
M.
,
Xiao
,
S.
,
Yang
,
Z.
, and
Sheng
,
P.
,
2014
, “
Acoustic Metasurface With Hybrid Resonances
,”
Nat. Mater.
,
13
(
9
), pp.
873
878
. 10.1038/nmat3994
35.
Zhu
,
H.
, and
Semperlotti
,
F.
,
2016
, “
Anomalous Refraction of Acoustic Guided Waves in Solids With Geometrically Tapered Metasurfaces
,”
Phys. Rev. Lett.
,
117
(
3
), p.
034302
. 10.1103/PhysRevLett.117.034302
36.
Su
,
X.
, and
Norris
,
A. N.
,
2016
, “
Focusing, Refraction, and Asymmetric Transmission of Elastic Waves in Solid Metamaterials With Aligned Parallel Gaps
,”
J. Acoust. Soc. Am.
,
139
(
6
), pp.
3386
3394
. 10.1121/1.4950770
37.
Shen
,
X.
,
Sun
,
C. T.
,
Barnhart
,
M. V.
, and
Huang
,
G.
,
2018
, “
Elastic Wave Manipulation by Using a Phase-Controlling Meta-Layer
,”
J. Appl. Phys.
,
123
(
9
), p.
091708
. 10.1063/1.4996018
38.
Su
,
X.
,
Lu
,
Z.
, and
Norris
,
A. N.
,
2018
, “
Elastic Metasurfaces for Splitting SV- and P-Waves in Elastic Solids
,”
J. Appl. Phys.
,
123
(
9
), p.
091701
. 10.1063/1.5007731
39.
Liu
,
Y.
,
Liang
,
Z.
,
Liu
,
F.
,
Diba
,
O.
,
Lamb
,
A.
, and
Li
,
J.
,
2017
, “
Source Illusion Devices for Flexural Lamb Waves Using Elastic Metasurfaces
,”
Phys. Rev. Lett.
,
119
(
3
), p.
034301
. 10.1103/PhysRevLett.119.034301
40.
Lee
,
H.
,
Lee
,
J. K.
,
Seung
,
H. M.
, and
Kim
,
Y. Y.
,
2018
, “
Mass-Stiffness Substructuring of an Elastic Metasurface for Full Transmission Beam Steering
,”
J. Mech. Phys. Solids
,
112
, pp.
577
593
. 10.1016/j.jmps.2017.11.025
41.
Cao
,
L.
,
Yang
,
Z.
,
Xu
,
Y.
, and
Assouar
,
B.
,
2018
, “
Deflecting Flexural Wave With High Transmission by Using Pillared Elastic Metasurface
,”
Smart Mater. Struct.
,
27
(
7
), p.
075051
. 10.1088/1361-665X/aaca51
42.
Cao
,
L.
,
Yang
,
Z.
,
Xu
,
Y.
,
Fan
,
S. W.
,
Zhu
,
Y.
,
Chen
,
Z.
,
Vincent
,
B.
, and
Assouar
,
B.
,
2020
, “
Disordered Elastic Metasurfaces
,”
Phys. Rev. Appl.
,
13
(
1
), p.
014054
. 10.1103/PhysRevApplied.13.014054
43.
Cao
,
L.
,
Yang
,
Z.
,
Xu
,
Y.
,
Chen
,
Z.
,
Zhu
,
Y.
,
Fan
,
S. W.
,
Donda
,
K.
,
Vincent
,
B.
, and
Assouar
,
B.
,
2021
, “
Pillared Elastic Metasurface With Constructive Interference for Flexural Wave Manipulation
,”
Mech. Syst. Sig. Process.
,
146
, p.
107035
. 10.1016/j.ymssp.2020.107035
44.
Lin
,
Z.
, and
Tol
,
S.
,
2020
, “
Elastic Metasurfaces for Low-Frequency Flexural Wavefront Control
,”
International Design Engineering Technical Conferences and Computers and Information in Engineering Conference
,
Virtual
,
Aug. 17–19
.
45.
Froehly
,
L.
,
Courvoisier
,
F.
,
Mathis
,
A.
,
Jacquot
,
M.
,
Furfaro
,
L.
,
Giust
,
R.
,
Lacourt
,
P. A.
, and
Dudley
,
J. M.
,
2011
, “
Arbitrary Accelerating Micron-Scale Caustic Beams in Two and Three Dimensions
,”
Opt. Exp.
,
19
(
17
), pp.
16455
16465
. 10.1364/OE.19.016455
46.
Zhang
,
P.
,
Li
,
T.
,
Zhu
,
J.
,
Zhu
,
X.
,
Yang
,
S.
,
Wang
,
Y.
,
Yin
,
X.
, and
Zhang
,
X.
,
2014
, “
Generation of Acoustic Self-Bending and Bottle Beams by Phase Engineering
,”
Nat. Commun.
,
5
(
1
), pp.
1
9
.
47.
Graff
,
K. F.
,
1991
, “Waves in Membranes, Thin Plates, and Shells,”
Wave Motion in Elastic Solids
,
Dover Publications, Inc.
,
New York
, pp.
213
272
.