Tunable piezoelectric metamaterial for Lamb waves using periodic shunted circuits

For this study, the $y$ direction width of the piezoelectric plate is assumed to be $w\gg\Lambda$, where $\Lambda$ is the width of one unit-cell, consisting of one open (uncovered) cell, and one electroded (covered) cell $\Lambda=a_{1}+a_{2}$. We assume that the thickness $h$ of the plate is significantly less than the width and the length of the entire plate $h\ll w$, $h\ll N\Lambda$, where $N$ is the number of unit-cells. We assume that the electrode thickness $T$ is much smaller than the plate thickness $T/h\ll 1$, as was the case in previous piezoelectric plate studies [9, 10].

We further narrow this study to certain piezoceramics with non-zero piezoelectric coefficients $e_{31}$, $e_{33}$, $e_{32}$, $e_{24}$, $e_{15}$, and permittivity $\varepsilon^{\mathrm{S}}_{11}=\varepsilon^{\mathrm{S}}_{22}\neq\varepsilon^{\mathrm{S}}_{33}$, which comes from poling in the thickness ($z$) direction. This is the case for piezoceramics commonly used for metamaterials [11, 9, 10] and transducers[13] (e.g. lead zirconium titanate (PZT) and barium titanate). We assume that the waves in the ($x,y$) plane are decoupled from the waves in the ($x,z$) plane [10], and only include waves in the ($x,z$) plane in our analysis.

For the electroded cells, which we herein call covered cells, we use a modified 2D version of the approximate model for a covered piezoelectric transducer detailed in [14, 15, 16]. This model was originally created for standalone transducers and composite piezoelectric/polymer transducers, and we extend it here for plates made entirely of piezoelectric material. We assume the dielectric displacement vector components $D_{1}=D_{2}=0$ and $\partial D_{3}/\partial x=\partial D_{3}/\partial z=0$. The dielectric displacement in the thickness $z$ direction is assumed to be $D_{3}=D_{0}e^{i\omega t}$. The displacement in the $x$ direction is assumed to be

$$\xi_{1}=\left[K_{1}\sin\left(\frac{\omega x}{v_{1}}\right)+K_{2}\cos\left(\frac{\omega x}{v_{1}}\right)\right]e^{i\omega t},$$

(10)

while displacement in the $z$ direction is assumed to be

$$\xi_{3}=\left[K_{3}\sin\left(\frac{\omega z}{v_{3}}\right)+K_{4}\cos\left(\frac{\omega z}{v_{3}}\right)\right]e^{i\omega t},$$

(11)

where $\omega$ is the angular frequency, and $K_{1-4}$ are constants. When equations (10) and (11) are used in the governing equations (1)–(3), along with the relation (5), the phase velocities are found to be

$\displaystyle v_{1}=\sqrt{\frac{c^{\mathrm{D}}_{11}}{\rho}};$

$\displaystyle v_{3}=\sqrt{\frac{c^{\mathrm{D}}_{33}}{\rho}},$

(12)

where

$\displaystyle c^{\mathrm{D}}_{11}=c^{\mathrm{E}}_{11}+\frac{e^{2}_{31}}{\varepsilon^{\mathrm{S}}_{33}};$

$\displaystyle c^{\mathrm{D}}_{33}=c^{\mathrm{E}}_{33}+\frac{e^{2}_{33}}{\varepsilon^{\mathrm{S}}_{33}}.$

(13)

A weak form of the free mechanical boundary condition is used to find constants $K_{1}$,$K_{2}$,$K_{3}$,$K_{4}$. The lengths $\Sigma_{1}$, $\Sigma_{2}$, $\Sigma_{3}$, and $\Sigma_{4}$ are along boundaries of a covered cell, as shown in Figure 1. The weak form integrals are

$$\int_{\Sigma_{1}}T_{11}(x=-a_{1}/2)\,d\Sigma_{1}=\int_{\Sigma_{2}}T_{11}(x=a_{1}/2)\,d\Sigma_{2}=0$$

(14)

$$\int_{\Sigma_{3}}T_{33}(z=-h/2)\,d\Sigma_{3}=\int_{\Sigma_{4}}T_{33}(z=h/2)\,d\Sigma_{4}=0.$$

(15)

Due to the lack of bending stress and the small width dimension $a_{1}$, terms with $e_{15}$ and $e_{24}$ are neglected. In order to find $\partial\phi/\partial x=0$, the (${x,z}$) plane piezoelectric coefficient is set to zero ($e_{31}=0$) [15]. This gives an expression for the electric field in the thickness direction [15]

$$E_{3}=\frac{D_{3}}{\varepsilon^{\mathrm{S}}_{33}}\left\{1-e_{33}\frac{\kappa\omega}{v_{1}}\left[\cos\left(\frac{\omega z}{v_{3}}\right)+\tan\left(\frac{\vartheta_{3}}{2}\right)\sin\left(\frac{\omega z}{v_{3}}\right)\right]\right\},$$

(16)

where

$$\kappa=\frac{v_{3}hc^{\mathrm{E}}_{11}h_{33}\omega a_{1}}{c^{\mathrm{E}}_{11}c^{\mathrm{D}}_{33}\omega^{2}ha_{1}-4(c^{\mathrm{E}}_{13})^{2}v_{1}v_{3}\tan(\vartheta_{1}/2)\tan(\vartheta_{3}/2)}.$$

(17)

The symbols and ratios $\vartheta_{1}=\omega a_{1}/v_{1}$, $\vartheta_{3}=\omega h/v_{3}$, and $h_{33}=e_{33}/\varepsilon^{\mathrm{S}}_{33}$ are kept consistent with [15]. The elastic velocity at the boundaries of the cell are

$\displaystyle\dot{\xi}_{1}(x=-a_{1})=u_{1};$

$\displaystyle\dot{\xi}_{1}(x=+a_{1})=u_{2}$

(18)

$\displaystyle\dot{\xi}_{3}(z=-h/2)=u_{3};$

$\displaystyle\dot{\xi}_{3}(z=+h/2)=u_{4}.$

(19)

External forces applied to the cell are related to the velocities using the following integration:

$\displaystyle\int_{\Sigma_{1}}T_{11}(x=-a_{1}/2)\,d\Sigma_{1}=-F_{1}$

$\displaystyle\int_{\Sigma_{2}}T_{11}(x=a_{1}/2)\,d\Sigma_{2}=-F_{2}$

(20)

$\displaystyle\int_{\Sigma_{3}}T_{33}(z=-h/2)\,d\Sigma_{3}=-F_{3}$

$\displaystyle\int_{\Sigma_{4}}T_{33}(z=h/2)\,d\Sigma_{4}=-F_{4}.$

(21)

This gives equations for the side forces per unit length

$$F_{1}=b_{1}u_{1}+b_{2}u_{2}+b_{5}(u_{3}+u_{4}),$$

(22)

$$F_{2}=b_{2}u_{1}+b_{1}u_{2}+b_{5}(u_{3}+u_{4}),$$

(23)

and top and bottom forces per unit length

$$F_{3}=b_{5}(u_{1}+u_{2})+b_{3}u_{3}+b_{4}u_{4},$$

(24)

$$F_{4}=b_{5}(u_{1}+u_{2})+b_{4}u_{3}+b_{3}u_{4}.$$

(25)

The voltage across the electrodes can be found by integrating (16) to give

$$V=b_{6}(u_{3}+u_{4})+Z_{\mathrm{eq}}I,$$

(26)

where $I$ is the current per unit length, and the electrical impedance of the plate is

$$Z_{\mathrm{eq}}=Z_{\mathrm{C0}}=\frac{1}{i\omega C_{0}},$$

(27)

with the clamped capacitance per unit length as

$$C_{0}=\frac{\varepsilon_{33}^{\mathrm{S}}a_{1}}{h}.$$

(28)

The coefficients in the above equations (22)–(26) are

$\displaystyle b_{1}=\frac{Z_{1}}{i\tan(\vartheta_{1})};$

$\displaystyle b_{2}=\frac{Z_{1}}{i\sin(\vartheta_{1})}$

(29)

$\displaystyle b_{3}=\frac{Z_{3}}{i\tan(\vartheta_{3})};$

$\displaystyle b_{4}=\frac{Z_{3}}{i\sin(\vartheta_{3})}$

(30)

$\displaystyle b_{5}=\frac{c_{13}^{\mathrm{E}}}{i\omega};$

$\displaystyle b_{6}=\frac{h_{33}}{i\omega},$

(31)

where the acoustic impedance per length in the $x$ direction is $Z_{1}=\rho hv_{1}$ and the acoustic impedance per length in the $z$ direction is $Z_{3}=\rho a_{1}v_{3}$. These equations differ from the equations formulated in [15] in that they are for 2D plane strain, and use forces per unit length. Equations (22)–(26) can be formulated into a $5\times 5$ impedance matrix corresponding to a 5 port element. When a shunted circuit with a complex impedance load $Z_{\mathrm{L}}$ is connected to the electrodes, The electrical impedance in equation (27) becomes

$$Z_{\mathrm{eq}}=\frac{Z_{\mathrm{L}}Z_{\mathrm{C0}}}{Z_{\mathrm{L}}+Z_{\mathrm{C0}}}.$$

(32)

In order for covered cells to couple with neighboring uncovered cells on both sides, the voltage port must be duplicated. As shown in the left-hand-side of Figure 2 (a), the 5 port covered cell element has only one voltage port, which is closed with load $Z^{j}_{\mathrm{L}}$. To create a cell with voltage ports on both sides, the one voltage port is duplicated, giving $V^{\mathrm{A}}_{1}$ and $V^{\mathrm{A}}_{2}$, which are the left side port voltage and the right side port voltage. The current per unit length from the left side is $I^{\mathrm{A}}_{1}$, while the current per unit length from the right side is $I^{\mathrm{A}}_{2}$. The voltages $V^{\mathrm{A}}_{1}$ and $V^{\mathrm{A}}_{2}$, along with forces $F_{i}$, $i=1,2,3,4,5,6$, constitute a 6 port element that can be connected on both the left and right side in a 1D network.

A 6 port covered cell element can be further simplified to a 4 port element when there is no material on the top or bottom of the plate that requires element modeling. The force at the bottom of the cell can be related to the velocity at the bottom of the cell by $F_{3}=Z_{\mathrm{B}}u_{3}$, where $Z_{\mathrm{B}}$ is the mechanical impedance of the layer below the plate. Likewise, the force at the top of the cell can be related to the velocity at the top of the cell by $F_{4}=Z_{\mathrm{T}}u_{4}$, where $Z_{\mathrm{T}}$ is the mechanical impedance of the layer above the plate. This eliminates the need for 2 of the 6 ports, giving the 4 port element shown in Figure 2(a). Equations (22)–(26) can be expressed in the impedance matrix form $q=Ar$, where $A$ is the impedance matrix, and $r$ is the velocity/current vector, and $q$ is the force/voltage vector. This impedance matrix formulation is fully expressed as

$$\begin{pmatrix}F_{1}\\ V_{1}\\ F_{2}\\ V_{2}\\ \end{pmatrix}=\begin{pmatrix}b_{1}-b_{5}^{2}b_{7}&b_{5}b_{6}b_{7}&b_{2}-b_{5}^{2}b_{7}&b_{5}b_{6}b_{7}\\ b_{5}b_{6}b_{7}&Z_{\mathrm{eq}}+b_{6}^{2}b_{7}&b_{5}b_{6}b_{7}&Z_{\mathrm{eq}}+b_{6}^{2}b_{7}\\ b_{2}-b_{5}^{2}b_{7}&b_{5}b_{6}b_{7}&b_{1}-b_{5}^{2}b_{7}&b_{5}b_{6}b_{7}\\ b_{5}b_{6}b_{7}&Z_{\mathrm{eq}}+b_{6}^{2}b_{7}&b_{5}b_{6}b_{7}&Z_{\mathrm{eq}}+b_{6}^{2}b_{7}\\ \end{pmatrix}\begin{pmatrix}u_{1}\\ I_{1}\\ u_{2}\\ I_{2}\\ \end{pmatrix}$$

(33)

with the use of an additional coefficient

$$b_{7}=\frac{2(b_{3}-b_{4})+Z_{\mathrm{B}}-Z_{\mathrm{T}}}{b_{3}^{2}-b_{4}^{2}+b_{4}Z_{\mathrm{B}}-b_{3}Z_{\mathrm{T}}}.$$

(34)

To formulate an impedance matrix equation for open (uncovered by electrodes) cells, electrostatic lumped impedances are implemented. Within open cells, there is a non-negligible component of the electric field in the $x$ direction [9]. We therefore assume a non-zero dielectric displacement in the $x$ direction $D_{1}$, and $\partial D_{1}/\partial x=\partial D_{3}/\partial z=0$ to satisfy (6). However, even if the effects of $D_{3}$ are neglected, equations like (22)–(26) cannot be formulated with the same assumed solutions for the displacements. Therefore, we simplify the problem by neglecting electro-mechanical coupling. This gives a linear impedance matrix

$$\mathbf{B}=\begin{pmatrix}b_{1}-b_{5}^{2}b_{7}&0&b_{2}-b_{5}^{2}b_{7}&0\\ 0&Z^{\mathrm{E}}_{11}&0&Z^{\mathrm{E}}_{12}\\ b_{2}-b_{5}^{2}b_{7}&0&b_{1}-b_{5}^{2}b_{7}&0\\ 0&Z^{\mathrm{E}}_{21}&0&Z^{\mathrm{E}}_{22}\\ \end{pmatrix}.$$

(35)

In the above matrix, null values $B_{12}=B_{14}=B_{21}=B_{23}=B_{32}=B_{34}=B_{41}=B_{43}=0$ are the electro-mechanical coupling terms. The potential difference between the top left corner and bottom left corner of the cell is $V^{\mathrm{B}}_{1}=Z^{\mathrm{E}}_{11}I^{\mathrm{B}}_{1}+Z^{\mathrm{E}}_{12}I^{\mathrm{B}}_{2}$, while the potential difference between the top right corner and bottom right corner of the cell is $V^{\mathrm{B}}_{2}=Z^{\mathrm{E}}_{21}I^{\mathrm{B}}_{1}+Z^{\mathrm{E}}_{22}I^{\mathrm{B}}_{2}$. The current per unit length going into the cell from the left is $I^{\mathrm{B}}_{1}$, while the current per unit length going into the cell from the right is $I^{\mathrm{B}}_{2}$. We calculate the electrical impedances $Z^{\mathrm{E}}_{11}$, $Z^{\mathrm{E}}_{12}$, $Z^{\mathrm{E}}_{21}$, and $Z^{\mathrm{E}}_{22}$ by estimating mutual capacitance per unit length.

The placement of mutual capacitance captures the first two modes of electro-mechanical propagation. The first mode of propagation has electric fields between parallel electrodes opposite in polarity with neighboring electrode pairs, as shown in Figure 3 (a). This occurs when $\lambda=\Lambda$ [9]. The capacitance between nearest neighbor top electrodes is labeled in Figure 3 (a) as $C_{1}$. To estimate $C_{1}$, we use an analytical method based on conformal mapping [17] which gives

$$C_{1}=\frac{\varepsilon_{\mathrm{ave}}}{4}\mathcal{M}(k_{\mathrm{T}}(\mu))$$

(36)

where $\varepsilon_{\mathrm{ave}}=(\varepsilon_{11}+\varepsilon_{33})/2$, and $\mu=a_{1}/a_{2}$. The function $\mathcal{M}$ is

$$\mathcal{M}\left(k(z,y)\right)\cong\begin{cases}\frac{2\pi}{\ln\left[2\frac{1+\left(1-k_{\mathrm{BCP}}^{2}\right)^{1/4}}{1-\left(1-k_{\mathrm{BCP}}^{2}\right)^{1/4}}\right]}&0\leq k_{\mathrm{BCP}}\leq\frac{1}{\sqrt{2}}\\ \frac{2}{\pi}\ln\left[2\frac{1+\sqrt{k_{\mathrm{BCP}}}}{1-\sqrt{k_{\mathrm{BCP}}}}\right]&\frac{1}{\sqrt{2}}\leq k_{\mathrm{BCP}}\leq 1\end{cases},$$

(37)

and $k_{\mathrm{S}}$ is a dimensionless geometric parameter calculated by

$$k_{\mathrm{S}}(z,y)=\frac{\sqrt{\sinh\left[\frac{\pi}{2}(z+y)\right]\sinh\left(\frac{\pi}{2}z\right)}}{\cosh\left[\frac{\pi}{2}\left(z+\frac{y}{2}\right)\right]}.$$

(38)

The capacitance between the top electrodes and the aligned bottom electrodes is partially taken into consideration by Eq. (28). However, this equation underestimates the capacitance between parallel plates in a continuous medium by neglecting fringing effects. To more accurately estimate parallel plate capacitance we use [18]

$$C_{\mathrm{p}}=\frac{1.15\varepsilon_{\mathrm{ave}}a_{1}}{h}+1.4\varepsilon_{\mathrm{ave}}(a_{1}+1)\left(\frac{2T}{h}\right)^{0.222}+1.03\varepsilon_{\mathrm{ave}}h\left(\frac{2T}{h}\right)^{0.728}$$

(39)

In the above equation, the first term is the parallel plate capacitance with a multiplicative factor of $1.15$, while the second and third terms account for fringe effects on both sides of the electrodes. The fringe capacitance $C_{2}$ within an open cell is taken as $C_{2}=(C_{\mathrm{p}}-C_{0})/2$.

The second mode of electromechanical propagation is when neighboring pairs of electrodes have electric fields with the same polarity, as shown in Figure 3 (b). This occurs when $\lambda\gg\Lambda$ [9]. In this mode, $C_{1}$ can be neglected. However, there is a non-negligable capacitance between top electrodes and the nearest neighbor bottom electrodes, labeled as $C_{3}$ in Figure 3 (b). We use a modified version of the parallel plate capacitance to account for the misalignment of the electrodes, which is

$$C_{3}=\frac{\varepsilon_{\mathrm{ave}}a_{1}}{d}+1.4\varepsilon_{\mathrm{ave}}(a_{1}+1)\left(\frac{2T}{d}\right)^{0.222}+1.03\varepsilon_{\mathrm{ave}}h\left(\frac{2T}{d}\right)^{0.728}$$

(40)

where $d=\sqrt{h^{2}+(a_{2}+a_{1})^{2}}$ is the adjusted distance between the electrodes.

Lumped impedances with both $C_{1}$ and $C_{3}$ accounts for both of these modes. The result is a lattice circuit, between two shunted impedance circuits, as shown in Figure 3 (c). The electrical impedance parameters are

$$Z^{\mathrm{E}}_{11}=\frac{Z_{3}Z_{2}^{2}(Z_{1}+2Z_{3})(Z1+Z_{3})+Z_{2}Z_{3}^{2}(Z_{1}^{2}+2Z_{1}Z_{3})}{2Z_{2}Z_{3}(Z_{1}+Z_{3})(Z_{1}+2Z_{3})+Z_{2}^{2}(Z_{1}+2Z_{3})^{2}+Z_{3}^{2}(Z_{1}^{2}+2Z_{1}Z_{3})};$$

(41)

$$Z^{\mathrm{E}}_{12}=\frac{Z_{2}^{2}(Z_{1}+2Z_{3})^{2}+Z_{3}^{2}(Z_{1}+2Z_{3})}{2Z_{2}Z_{3}(Z_{1}+2Z_{3})(Z_{1}+Z_{3})+Z_{2}^{2}(Z_{1}+2Z_{3})^{2}+Z_{3}^{2}(Z_{1}^{2}+2Z_{1}Z_{3})};$$

(42)

where $Z_{1}=1/(i\omega C_{1})$, $Z_{2}=1/(i\omega C_{2})$, and $Z_{3}=1/(i\omega C_{3})$. The other terms in the symmetric matrix are $Z^{\mathrm{E}}_{22}=Z^{\mathrm{E}}_{11}$ and $Z^{\mathrm{E}}_{21}=Z^{\mathrm{E}}_{12}$.

To calculate scattering parameters for plates consisting of multiple electrode strips, scattering matrices are used. Scattering matrices are commonly used for the transfer matrix method when modeling electromagnetic wave propagation, due to being numerically stable and intrinsically containing the refection and transmission coefficients [19]. One can transform the impedance matrix of a covered cell $\mathbf{A}$ into a scattering matrix $\mathbf{S}^{\mathrm{A}}$ with the formula [20]

$$\mathbf{S}^{\mathrm{A}}=\mathbf{G}_{\mathrm{ref}}(\mathbf{A}-\mathbf{Z}^{*}_{\mathrm{ref}})(\mathbf{A}+\mathbf{Z}_{\mathrm{ref}})^{-1}\mathbf{G}^{-1}_{\mathrm{ref}}.$$

(43)

The reference impedance matrices are $Z_{\mathrm{ref},ij}=Z_{0,ij}\delta_{ij}$, and $G_{\mathrm{ref},ij}=1/\sqrt{|Z_{0,ij}|}\delta_{ij}$, where $\delta_{ij}$ is the Kronecker delta function. The mechanical reference impedances are assumed to be $Z_{0,11}=Z_{0,33}=\rho hv_{1}$, while the electrical reference impedances are assumed to be $Z_{0,22}=Z_{0,44}=50$ $\Omega$. Likewise, one can transform the impedance matrix of an uncovered cell $\mathbf{B}$ into a scattering matrix $\mathbf{S}^{\mathrm{B}}$.

The S matrix for a symmetric unit cell consists of one uncovered cell and one covered cell, shown in Figure 2 (c). The scattering matrix of a unit cell can be found with

$$\mathbf{S}^{\mathrm{UC}}=\mathbf{S}^{\mathrm{A}}\otimes\mathbf{S}^{\mathrm{B}},$$

(44)

where $\otimes$ is the Redheffer star product [19], which is detailed in Appendix A. The scattering matrix $\mathbf{S}^{\mathrm{UC},j}$ for the $j$th unit cell gives the relation

$$\begin{pmatrix}F^{j,-}_{1}\\ V^{j,-}_{1}\\ F^{j,+}_{2}\\ V^{j,+}_{2}\\ \end{pmatrix}=\begin{pmatrix}S^{j}_{11}&S^{j}_{12}&S^{j}_{13}&S^{j}_{14}\\ S^{j}_{21}&S^{j}_{22}&S^{j}_{23}&S^{j}_{24}\\ S^{j}_{31}&S^{j}_{32}&S^{j}_{33}&S^{j}_{34}\\ S^{j}_{41}&S^{j}_{42}&S^{j}_{43}&S^{j}_{44}\\ \end{pmatrix}\begin{pmatrix}F^{j,+}_{1}\\ V^{j,+}_{1}\\ F^{j,-}_{2}\\ V^{j,-}_{2}\\ \end{pmatrix}$$

(45)

where superscript $+$ and $-$ denote frontward and backward waves, respectively.

As shown in Figure 2 (c), a scattering matrix for a section of plate with $N$ electrode pairs can be found by coupling $N-1$ unit cells with a covered cell. This global scattering matrix can be found with

$$\mathbf{S}^{\mathrm{G}}=\mathbf{S}^{\mathrm{UC},1}\otimes\mathbf{S}^{\mathrm{UC},2}\otimes\dots\mathbf{S}^{\mathrm{UC},N-1}\otimes\mathbf{S}^{\mathrm{A}}.$$

(46)

If the plate has reflection and transmission regions without electrodes, asymmetric scattering matrices $\mathbf{S}^{\mathrm{RFL}}$ and $\mathbf{S}^{\mathrm{TRS}}$ coupling the cells to the reflection and transmission regions can be used [19]. Cells in this 2D plate model can each have different electrical and mechanical boundary conditions, as well as different material properties and thickness.

The dispersion relations for infinitely periodic plates can also be evaluated with the scattering matrix of one unit cell $j$. Terms in the scattering matrix $\mathbf{S}^{\mathrm{UC}}$ can be moved around to form a relation between wave amplitudes from unit-cell ports $F^{j,+}_{2}$, $V^{j,+}_{2}$ and ports $F^{j,+}_{1}$, $V^{j,+}_{1}$, given as

$$\begin{pmatrix}0&0&-S^{\mathrm{UC}}_{13}&-S^{\mathrm{UC}}_{14}\\ 0&0&-S^{\mathrm{UC}}_{23}&-S^{\mathrm{UC}}_{24}\\ 1&0&-S^{\mathrm{UC}}_{33}&-S^{\mathrm{UC}}_{34}\\ 0&1&-S^{\mathrm{UC}}_{43}&-S^{\mathrm{UC}}_{44}\\ \end{pmatrix}\begin{pmatrix}F^{j,+}_{2}\\ V^{j,+}_{2}\\ F^{j,-}_{2}\\ V^{j,-}_{2}\\ \end{pmatrix}=\begin{pmatrix}S^{\mathrm{UC}}_{11}&S^{\mathrm{UC}}_{12}&-1&0\\ S^{\mathrm{UC}}_{21}&S^{\mathrm{UC}}_{22}&0&-1\\ S^{\mathrm{UC}}_{31}&S^{\mathrm{UC}}_{32}&0&0\\ S^{\mathrm{UC}}_{41}&S^{\mathrm{UC}}_{42}&0&0\\ \end{pmatrix}\begin{pmatrix}F^{j,+}_{1}\\ V^{j,+}_{1}\\ F^{j,-}_{1}\\ V^{j,-}_{1}\\ \end{pmatrix}.$$

(47)

The Bloch wave theory for infinitely periodic plates then gives the relation [21]

$$\begin{pmatrix}F^{j,+}_{2}\\ V^{j,+}_{2}\\ F^{j,-}_{2}\\ V^{j,-}_{2}\\ \end{pmatrix}=e^{-i\beta\Lambda}\begin{pmatrix}F^{j,+}_{1}\\ V^{j,+}_{1}\\ F^{j,-}_{1}\\ V^{j,-}_{1}\\ \end{pmatrix}$$

(48)

where $\beta$ is the effective wavenumber at a given frequency. This relation allows us to formulate

$$\begin{pmatrix}0&0&-S^{\mathrm{UC}}_{13}&-S^{\mathrm{UC}}_{14}\\ 0&0&-S^{\mathrm{UC}}_{23}&-S^{\mathrm{UC}}_{24}\\ 1&0&-S^{\mathrm{UC}}_{33}&-S^{\mathrm{UC}}_{34}\\ 0&1&-S^{\mathrm{UC}}_{43}&-S^{\mathrm{UC}}_{44}\\ \end{pmatrix}\begin{pmatrix}F^{j,+}_{1}\\ V^{j,+}_{1}\\ F^{j,-}_{1}\\ V^{j,-}_{1}\\ \end{pmatrix}=e^{-i\beta\Lambda}\begin{pmatrix}S^{\mathrm{UC}}_{11}&S^{\mathrm{UC}}_{12}&-1&0\\ S^{\mathrm{UC}}_{21}&S^{\mathrm{UC}}_{22}&0&-1\\ S^{\mathrm{UC}}_{31}&S^{\mathrm{UC}}_{32}&0&0\\ S^{\mathrm{UC}}_{41}&S^{\mathrm{UC}}_{42}&0&0\\ \end{pmatrix}\begin{pmatrix}F^{j,+}_{1}\\ V^{j,+}_{1}\\ F^{j,-}_{1}\\ V^{j,-}_{1}\\ \end{pmatrix}$$

(49)

which is in the generalized eigenvalue form. This eigenvalue problem can be evaluated to find $\lambda=e^{-i\beta\Lambda}$ for a given frequency. The effective wavenumber $\beta$ can then be found with

$$\cos(\beta\Lambda)=\frac{1}{2}\left(e^{-i\beta\Lambda}+e^{i\beta\Lambda}\right)$$

(50)

Given the scattering parameters for a unit-cell or for an entire plate consisting of multiple cells, the effective elastic impedance and effective elastic wavenumber can be calculated. We use a method of extracting the effective mechanical impedance that is commonly used for electromagnetic metamaterials [22]. The effective mechanical impedance can be found from the mechanical scattering coefficients with [22]

$$Z_{\mathrm{eff}}=\pm Z_{0}\sqrt{\frac{(S_{11}+1)^{2}-S_{31}^{2}}{(S_{11}-1)^{2}-S_{31}^{2}}}$$

(51)

where $Z_{0}$ is the reference impedance, which we set as equal to the mechanical impedance per length $Z_{0}=Z_{0,11}$.

The effective wavenumber $k_{\mathrm{eff}}$ is found by the definition $P_{\mathrm{A}}=e^{ik_{\mathrm{eff}}N\Lambda}$. The complex value $P_{\mathrm{A}}$ can be found at a given frequency by [22]

$$P_{\mathrm{A}}=\frac{S_{31}(Z_{\mathrm{eff}}+Z_{0})}{(Z_{\mathrm{eff}}+Z_{0})-S_{11}(Z_{\mathrm{eff}}-Z_{0})}$$

(52)

The effective wavenumber is then calculated by [22]

$$k_{\mathrm{eff}}=\frac{1}{N\Lambda}\left[-\left(\arg(P_{\mathrm{A}})+2\pi q\right)+i\log|P_{\mathrm{A}}|\right];$$

(53)

where $q$ is the branch index that must be chosen for each frequency.

Eigenfrequency FEA studies were conducted with and without altered property values to estimate the effect of some of the key assumptions made in the analytical model formulation. Studies with and without certain piezoelectric coefficients set to zero ($e_{15}=e_{31}=0$) were conducted. A 2D rectangular piezoelectric solid with perfectly conductive line electrodes on the top and bottom was used. A triangular mesh with a maximum dimension $\delta=\lambda_{\mathrm{m}}/200$, where $\lambda_{\mathrm{m}}$ is the smallest wavelength in the study, was utilized to divide the occupied geometry. Shunted electrical circuits with passive elements were accounted for by lumped parameter differential equations coupled to the partial differential equations for the piezoelectric solid. The Bloch wave condition was used for the right and left boundaries of a unit-cell, allowing the calculation of eigenfrequencies for an infinitely periodic plate. A pair of band diagrams for $h=2$ mm, $a_{1}=1$ mm, $a_{2}=2$ mm, and $L=0.1$ mH, with and without $e_{15}=e_{31}=0$ is shown in Figure 4.

As shown in Figure 4, the largest impact in the dispersion relation made by setting $e_{31}=e_{15}=0$ is found in the zeroth order asymmetric lamb waves $A_{0}$. This is due to the absence of the $\epsilon_{31}$ and $\epsilon_{15}$ strain terms. The analysis with the neglected terms estimates a higher low-frequency hybridization band gap for $S_{0}$, with a gap that spans 151–177 kHz, compared to 137–169 kHz without any terms neglected. Neglecting piezoelectric coefficients ($e_{31}=e_{15}=0$) produces a slight change in the hybridization band gap, but does not change the overall shape and structure of the symmetric mode $S_{0}$ band diagram.

Methods in Section III were used to numerically calculate band diagrams for infinitely periodic plates. It was assumed that the plates were mechanically free on the top and bottom ($Z_{\mathrm{T}}=0$, $Z_{\mathrm{B}}=0$). These band diagrams were compared with band diagrams calculated using full FEA with $e_{15}\neq 0$ and $e_{31}\neq 0$. Line electrode elements were used in the mesh for the FEA calculations. A periodic condition was used for both the mechanical and electrostatic boundary on the left and right-hand side of the geometry. Equations (49) and (50) were solved numerically, giving effective elastic wavenumbers $\beta$ at frequencies $f\in{0.01,0.7}$ MHz, where $f=\omega/(2\pi)$. These solutions are shown with red $+$ markers in Figure 5 for $h=2$ mm and $a_{1}=1$ mm, $a_{2}=2$ mm, and L=$0.1$ mH, while the FEA solutions are shown with black dots.

As shown in Figure 5, the analytical results compare better with the FEA results near the hybridization band gap than the high frequency Bloch wave band gap. The analytical method estimates the Bloch wave band gap to be 620–625 kHz, while the FEA results estimate it to be 568–574 kHz. We define the percent difference in the middle frequency of the hybridization bandgap as

$$\Delta f_{\mathrm{m}}=100\%\times\frac{f^{\mathrm{A}}_{\mathrm{m}}-f^{\mathrm{F}}_{\mathrm{m}}}{f^{\mathrm{F}}_{\mathrm{m}}},$$

(54)

where subscripts $A$ and $F$ refer to the analytical method and the FEA results, respectively. For this example, $\Delta f_{\mathrm{m}}=-6.7$kHz (-4.4%). The percent difference in the width of the hybridization band gap $\Delta f_{\mathrm{g}}$, which is calculated with Eq. (54) using $f^{\mathrm{A}}_{\mathrm{g}}$ and $f^{\mathrm{F}}_{\mathrm{g}}$, in this example is 11.9 kHz (-37%). The eigenfrequencies at $k=1$ ($f_{\mathrm{l}}$) and $k=0$ ($f_{\mathrm{u}}$) also show some disagreement. The percent differences $\Delta f_{\mathrm{u}}$ and $\Delta f_{\mathrm{l}}$ were calculated using Eq. (54) with $f^{\mathrm{A}}_{\mathrm{u}}$,$f^{\mathrm{F}}_{\mathrm{u}}$,$f^{\mathrm{A}}_{\mathrm{l}}$, and $f^{\mathrm{F}}_{\mathrm{l}}$. For the example geometry used for Figure 5, $\Delta f_{\mathrm{l}}=2.7\%$ and $\Delta f_{\mathrm{u}}=-7.2\%$. The analytical method underestimates the band gap by placing it at lower frequencies than those found from FEA, and underestimates the width of the band gap, but does very well capturing the electromechanical coupling, the overall band structure, and the hybridization band gap.

The example band diagrams in Figure 4 and Figure 5, show that while the disagreement between the analytical method and FEA is partially due to the neglected piezoelectric coefficients $e_{31}=e_{15}=0$, it is also due to other factors. The FEA band diagram in 4 is at slightly higher freqencies when $e_{31}=e_{15}=0$. However, the band diagram calculated using the approximate analytical method is at slightly lower frequencies than the band diagram from FEA. The decoupling of the electrical and mechanical impedances in Equation (35), and the use of an approximate electrical capacitance for electrical impedance (Eq. (36)–(40) also seem to cause discrepancies between band diagrams calculated using FEA and the analytical method.

In order to understand the efficacy of the analytical model in predicting what frequencies the hybridization band gap occurs, the percent difference $\Delta f_{\mathrm{m}}$ was calculated for different plate thickness values $h$ and different gap values between the electrodes $a_{2}$. In Figure 6 (a), the $\Delta f_{\mathrm{m}}$ values are shown as percentages for $a_{2}\in\{0.8,...,3.2\}$mm and $h\in\{2,3,4\}$mm. In Figure 6 (b)-(d), It can be seen from Figure 6 (a) that $\Delta f_{\mathrm{m}}$ is greater for larger $h$ values, and less for lower $a_{2}$ values. The $|\Delta f_{\mathrm{m}}|$ values are all below $13\%$ for these geometries, and are as low as $0.45\%$, demonstrating a close agreement between the approximate analytical method and multiphysics FEA for these geometries. The fact that $\Delta f_{\mathrm{m}}$ is not consistently positive likely has to due with the accuracy of the capacitance estimations and/or the use of Eq. (41 and (42) for estimating $Z^{\mathrm{E}}_{11}$ and $Z^{\mathrm{E}}_{12}$.

While the analytical method closely agrees with the FEA calculated $f_{\mathrm{m}}$, the band gap width $f_{\mathrm{g}}$, and the eigenfrequencies $f_{\mathrm{u}}$, $f_{\mathrm{l}}$ show greater disagreement. The analytical method underestimates $f_{\mathrm{u}}$ by $5$-$20\%$. The analytical method overestimates the eigenfrequency of the $S_{0}$ mode at $k=1$ ($f_{\mathrm{l}}$) by $0.1$-$47\%$. For lower values of $h$, and higher values of $a_{2}$, $|\Delta f_{\mathrm{l}}|$ and $|\Delta f_{\mathrm{u}}|$ are reduced. Figure 6 (b) shows how the approximate analytical method underestimates the width of the band gap by $19$-$56\%$. This underestimation also reduces with an increase in $a_{2}$, but reduces with a larger $h$. These calculations demonstrate that the agreement between the analytical method and the FEA method depends greatly on the thickness $h$ and electrode gap $a_{2}$.

To investigate the effect of the electrical capacitance equations on calculated band diagrams, calculated capacitances from FEA steady state electrostatic simulations were compared with capacitances calculated using equations (36)–(40). Figure 7 (a)–(c) shows plots of the capacitance values for a set of width values $a_{2}\in\{0.2,...,4.0\}$.

Figure (7) shows a close agreement in capacitance values between the analytical methods and the FEA simulation. In Figure (7) (a) it can be seen that the capacitance between parallel plates for this geometry, calculated using Eq. (39) to include fringing effects is nearly double the capacitance calculated using the simple parallel plate equation (28). The capacitance $C_{\mathrm{p}}$ is within 3.4% of the capacitance values from the simulation for all values of $a_{2}$. The top electrode nearest neighbor capacitance $C_{1}$ is within 13% of the capacitance values from the simulation for all values of $a_{2}$. The greater discrepency in $C_{1}$ for $a_{2}>2$mm is due to the limitation of the method for calculated capacitance when the distance between electrodes exceeds the thickness of the material. The capacitance between top and nearest neighbor bottom electrodes $C_{3}$ agreed well with simulated values, staying within 5.7% of the FEA values for all values of $a_{2}$. Based on how well the FEA calculated capacitances compare with the analytically calculated capacitances, it is likely that innacuracies in computed capacitances contribute to some, but not all, of the band diagram discrepencies between analytical and FEA methods.

Finite plate models with $N=6$ electrodes were used to compute global scattering coefficients with Equation (46). The example geometry of $h=2$ mm, $a_{1}=0.5$ mm, $a_{2}=0.5$ mm was used with $L=1$ mH to avoid Fabry Perot resonances in the plate and to ensure $\lambda\ll\Lambda$ near the hybridization band gap. The reflection coefficient $S^{\mathrm{G}}_{11}$ and the transmission coefficient $S^{\mathrm{G}}_{13}$ for the acoustic wave is shown for frequencies $f\in\{0.01,0.1\}$ MHz in Figure 8 (a). The reflection coefficient $S^{\mathrm{G}}_{22}$ and the transmission coefficient $S^{\mathrm{G}}_{24}$ for the electric potential wave is shown in Figure 8 (b) below. While the computed global reflection coefficients change with the number of unit-cells $N$, Figure 8 gives an example of what the scattering coefficient look like for band gaps near $\lambda\ll\Lambda$ with a practical amount of electrode pairs in the array.

The scattering coefficients in Figure 8 have multiple inflection frequencies, associated with different propagation modes. Near 42 kHz, there is a decrease in $S^{\mathrm{G}}_{24}$ and $S^{\mathrm{G}}_{13}$. This is associated with the first mode of propagation, shown symbolically by the mutual capacitances in Figure 3 (a). The decrease in transmission at 57 kHz is associated with the second mode of propagation, which is shown symbolically by the mutual capacitances in Figure 3 (b). For this geometry, this occurs at a slightly higher frequency than the first mode, and is more pronounced for a $S_{0}$ mode excitation.

To show how the scattering parameters change with varying the inductance $L$ and by varying the height $h$, the middle of the second mode inflection corresponding to a band gap $f_{\mathrm{m}}$ was calculated, as well as the quality factor $Q$. The quality factor is defined as $Q=f_{\mathrm{m}}/(f_{\mathrm{w1}}-f_{\mathrm{w2}})$, where $f_{\mathrm{w1}}$ and $f_{\mathrm{w2}}$ are the first and second frequencies at half of the decrease in transmission $S^{\mathrm{G}}_{13}$. Figure 10 (a) shows $f_{\mathrm{m}}$ and $Q$ for $a_{1}=0.5$ mm, $a_{2}=0.5$ mm, $L=1$ mH, and $h\in\{2$,$3$,$4\}$ mm, while Figure 10 (b) shows $f_{\mathrm{m}}$ and $Q$ for $a_{1}=0.5$ mm, $a_{2}=0.5$ mm, $h=2$ mm, and $L\in\{1,10,100,1000\}$ mH.

The band gap moves to higher frequencies with an increase in $h$. This is shown in Figure 9 (a), with $f_{\mathrm{m}}=56.6$kHz for $h=2$mm, and $f_{\mathrm{m}}=73.5$kHz for $h=4$mm. The quality factor $Q$, however, gets smaller with an increase in $h$. This is to be expected, due to the electrical boundary conditions having a less pronounced effect on the scattering parameters as $h$ is increased.

The band gap moves to lower frequencies with an increase in inductance $L$. This is shown in Figure 9 (b), with $f_{\mathrm{m}}=180$kHz for $L=10^{-4}$ H, and $f_{\mathrm{m}}=5.65$kHz for $L=10^{-1}$mm. The quality factor $Q$ also increases with an increase in $L$. Therefore, if high $Q$ is to be desired for an application, increasing $L$ until the desired band gap frequency is reached would be the simplest design consideration.

The frequency shift in $f_{\mathrm{m}}$ due to inductance is expected and predicted by an infinite circuit model of this problem. If we consider the infinite plate model, the hybridization band gap is an avoided crossing with the purely electrical circuit band and the $S_{0}$ mode lamb wave band. A simplified electrical circuit for the first electrical mode (Figure 3 (a)) used in [9, 10] is

$$\sin(\Lambda/2)=\sqrt{\frac{1}{LC_{1}\omega^{2}}-\frac{C_{\mathrm{p}}}{C_{1}}}$$

(55)

From this equation, it can be predicted that the electrical band will shift to lower frequencies with an increase in $L$ or $C_{1}$. The upper frequency of the band gap is predicted in the above equation by the LC resonance, which for $h=2$ mm, $a_{1}=1$ mm, and $a_{2}=2$ mm is $1/\sqrt{LC_{0}}=203$ kHz. The lower frequency is predicted to be $95.4$ kHz. The approximate analytical method presented in this article gives $f^{\mathrm{A}}_{\mathrm{u}}=157$ kHz and $f^{\mathrm{A}}_{\mathrm{l}}=122$ kHz, and the FEA gives $f^{\mathrm{F}}_{\mathrm{u}}=170$ kHz and $f^{\mathrm{F}}_{\mathrm{l}}=119$ kHz. Though the above equation (55) gives larger differences in band gap estimation with FEA than the analyical method presented in this article, it may be useful for building intuition on the affect of $L$ and geometry that changes $C_{1}$ and $C_{0}$, and for estimating the middle of the band gap $f_{\mathrm{m}}$.