Mode couplings associated with elastic wave propagation through three-dimensional multiplex structures, as manifested by asymmetric eigenmodes and dissipation, determine the efficiency of electromechanical structures. As a result, it is critical to predict electroelastic symmetric modes such as thickness expander and radial modes, as well as asymmetric flexural modes, while accounting for material losses. Multiplex electromechanical structures include multi-layered through-wall ultrasound power transfer (TWUPT) systems. Physical processes that support TWUPT include vibrations at a transmitting/acoustic source element, elastic wave propagation through a barrier and coupling layers, piezoelectric transduction of elastic vibrations at a receiving element, and spatial resonances of the transmitting and receiving elements. We investigate mode couplings in an optimized modal TWUPT system, including their physical origins, models used to describe them, and regimes of weak and strong couplings. The system layout optimization is defined in terms of size (volume), operating frequency, and matching circuit load optimization. A computational model is developed and utilized in conjunction with experimental modal characterization to highlight the impact of eigenmode features on optimization results. Several behavioral modes are identified and analyzed. The interaction of symmetric radial and asymmetric flexural modes causes the system damping to increase and the device's overall efficiency to decrease. The electromechanical coupling factor value is likewise reduced as a result of this. Such occurrences are explained by the flow of energy between modes as they interact. The present work also proposes design guidelines to improve the performance of TWUPT systems based on exploiting inherent physical phenomena.

## I. INTRODUCTION

The electroelastic mode coupling phenomena have gained the interest of researchers over the last few decades due to their intriguing and yet complicated dynamics.^{1} These phenomena are characterized by energy exchange between two vibrating modes in which energy transfers from one to the other. Due to their complexity, vibration mode couplings in multiplex electromechanical structures have not been fully explored. Among others, multi-layered through-wall ultrasound power transfer (TWUPT) systems^{2–6} constitute one type of these structures that could be used to study the mode couplings phenomena. The TWUPT system, as shown in Fig. 1, consists of vibrating a transmitting piezoelectric element, propagating elastic waves through an elastic wall and bonded layers, and piezoelectric transduction of elastic vibrations at a receiving disk element. TWUPT systems are a new potential technology for wireless power transmission to hard-to-reach components^{2} in critical engineering applications ranging from ultrasound nondestructive evaluation^{7} to medical ultrasound, where acoustic wave manipulation and focusing to confine acoustic intensity is of great importance.^{8–10}

Mode couplings are achieved through mechanical, electrical, or internal resonance within the structure.^{11} Such phenomenon is frequently observed in MEMS structures and has been exploited to enhance different functionalities. Hajjaj *et al*.^{11} demonstrated a bandpass filter with sharp roll off from the passband to the stopband by exploiting nonlinear softening, hardening, and veering phenomena that occur when the frequencies of two vibration modes approach each other. The results of their research show that veering occurs when both modes become connected and exchange energy as a result. They also demonstrated that quick roll off from the passband to the stopband can be achieved by driving both modes nonlinearly and electrostatically near the veering regime, resulting in softening and hardening behavior in the first and third modes, respectively. The method for stabilizing oscillation frequency through mode coupling by the internal resonance condition is provided by Antonio *et al*.^{12} Their proposed method can be applicable to a wide range of micro- and nanomechanical oscillators. They also found that the coupling can be achieved not only with the torsional mode but also with the out-of-plane flexural mode in other comparable resonators, making the design simple to apply in other resonators.

Internal resonances within the structure itself can be used to achieve coupling in electromechanical systems.^{11} Internal resonance occurs when the ratio of the frequencies associated with two interacting vibration modes of a system becomes an integer ratio such as 1:1,^{13} 1:2,^{14} 1:3,^{12} 2:1,^{15} and 3:1.^{13} Nonlinear systems exhibit mode couplings via internal resonances with varying natural frequency ratios. The presence of such an internal resonance ratio in nonlinear systems is influenced by their nonlinearity degrees. For example, in a system with a 2:1 internal resonance, the natural frequencies of two resonant modes caused by quadratic nonlinearities have a frequency ratio of 2:1. Nonlinearity allows a resonated mode's superharmonic or subharmonic frequency to coincide with another mode of the system. Ouakad *et al*.^{13} studied the 1:1 internal resonance between the first and second asymmetric modes of an arch-shaped microstructure with varying initial rise values. Zhang *et al*.^{15} investigated the 2:1 internal resonance in a micromechanical cantilever beam resonator between extensional and flexural modes. The effect of the structure size and aspect ratio of the electroded and unelectroded regions on mode coupling effects for partially electrode thin-film bulk acoustic wave resonators was investigated by Li *et al*.^{16} Linear TWUPT systems exhibit mode coupling effects via 1:1 internal resonance, which has a significant impact on their overall efficiency. As such, mode coupling through internal resonances in TWUPT systems will be investigated and described in this work.

For understanding the elements of mode couplings in TWUPT systems, we propose experimental modal characterization, a geometric-dependent framework, and a 3D optimization approach. The varying parameters of the optimization problem include the size (volume), operating frequency, and matching circuit load. The emphasis is on the use of the genetic algorithm (GA) for characterization as well as the gradient-based method in conjunction with a finite element model (FEM) to maximize objective functions such as efficiency and electrical output power. A variety of numerical examples are used to demonstrate the effects of eigenmode qualities on optimization results, which are then validated using experimental measurements. We may examine the feasibility of running systems with a wide range of transducer aspect ratios in a number of vibration modes using the proposed models and experiments. However, one-dimensional modeling approaches such as analytical and equivalent circuity modeling,^{2,17,18} which were previously proposed, could not be employed to undertake such a study. This work sheds light on design ideas in the context of TWUPT optimization. The paper concludes with a discussion of possible issues related to the operation of TWUPT systems and potential future research.

## II. ANALYSIS OF MODE COUPLINGS IN ELECTRO-MECHANICAL TWUPT SYSTEMS THROUGH GEOMETRIC OPTIMIZATION

### A. Modal characterization and parameter identification using electrical impedance measurement in conjunction with genetic algorithm

TWUPT systems, shown in Fig. 1, are made up of an elastic barrier surrounded by piezoelectric transmitting and receiving elements (PZT, lead zirconate titanate). Elastic wave propagation takes place in the barrier and coupling layers. Inverse and direct piezoelectric transduction of elastic vibrations occurs at the transmitting and receiving piezoelectric elements. An electrical resistive load is connected at the receiver end where the electric output power and efficiency are calculated.

Electrical impedance measurements are performed for a freely hanging disk (PZT-4/SM111 from STEMiNC) using an HP4192A impedance analyzer throughout a recurrence run of 0–600 kHz to characterize the piezoelectric transmitter and receiver disks in Fig. 1. A high-fidelity FEM computational model implemented in COMSOL is used to predict the system performance under different excitation and attached resistive loads. The model is communicated directly with GA optimization function via MATLAB LiveLink. The data are utilized to create the GA fitness function, which evaluates the difference between the prospective and experimental resonance (local minimum) and anti-resonance (local maximum) frequencies, i.e., short- and open-circuit frequencies, as well as their amplitudes. Table I lists the ten independent variables used in the process, namely, the elastic constants, piezoelectric constants, and dielectric coefficients. We note that multiple GA runs are undertaken with different population sizes to ensure that the global minimum point was solved. As a result, in Figs. 2(a) and 2(b), the measured impedance frequency response functions (FRFs) are aligned with FEM's predictions, and the modified material properties are listed in Table I. In addition, the frequency-varying material properties used in this study are described in the supplementary material (see Fig. S1). It should be noted that our tests, carried out at the modest excitation level (*V _{in}* ≈ 1 V), demonstrate the TWUPT system's linear properties, where the output voltage signal's amplitudes are proportional to the input voltage amplitudes. We note that the strains generated for 1 V voltage amplitude are not significant to trigger any material nonlinear behavior in the layers of the TWUPT system. As a result, the FEM was performed assuming linear piezoelectricity and neglecting elastic coupling and dissipative nonlinearities.

Piezoelectric disks | ||||||||||

Density (kg m^{−3})
. | Elastic constants (10^{10} Pa)
. | Piezoelectric constants (C m^{−2})
. | Dielectric coefficients . | |||||||
---|---|---|---|---|---|---|---|---|---|---|

$\rho $ | $c11E$ | $c12E$ | $c13E$ | $c33E$ | $c44E$ | e_{31} | e_{33} | e_{15} | $\u03f511S$ | $\u03f533S$ |

7650 | 15 | 7.8 | 8.05 | 13.1 | 2.564 | −4.5 | 14.4 | 12 | 680ε_{0} | 550ε_{0} |

Piezoelectric disks | ||||||||||

Density (kg m^{−3})
. | Elastic constants (10^{10} Pa)
. | Piezoelectric constants (C m^{−2})
. | Dielectric coefficients . | |||||||
---|---|---|---|---|---|---|---|---|---|---|

$\rho $ | $c11E$ | $c12E$ | $c13E$ | $c33E$ | $c44E$ | e_{31} | e_{33} | e_{15} | $\u03f511S$ | $\u03f533S$ |

7650 | 15 | 7.8 | 8.05 | 13.1 | 2.564 | −4.5 | 14.4 | 12 | 680ε_{0} | 550ε_{0} |

Aluminum wall . | ||||
---|---|---|---|---|

Density (kg m^{−3})
. | Elastic constants (10^{9} Pa)
. | Poisson's ratio . | Quality factor . | |

$\rho $ | E | v | Q_{m} | Q_{d} |

2710 | 69 | 0.33 | 1000 | 750 |

Aluminum wall . | ||||
---|---|---|---|---|

Density (kg m^{−3})
. | Elastic constants (10^{9} Pa)
. | Poisson's ratio . | Quality factor . | |

$\rho $ | E | v | Q_{m} | Q_{d} |

2710 | 69 | 0.33 | 1000 | 750 |

Bonding layer (epoxy) | |||

1150 | 4 (average) | 0.33 | 8.4 |

Bonding layer (epoxy) | |||

1150 | 4 (average) | 0.33 | 8.4 |

where *ε*_{0} = 8.854 × 10^{−12} farads m^{−1}. The elastic under constant electric field E(*c*^{E}), piezoelectric (*e*), and dielectric under constant strain S(*ε*^{S}) matrices are given in the Appendix.

In TWUPT systems, there are three types of energy losses to consider: dielectric, elastic, and piezoelectric losses. To account for these losses, Holland^{19} was the pioneer to employ an imaginary portion for the piezoelectric constants. Uchino *et al*.^{20} proposed later to divide each loss into intensive and extensive components. Several studies^{21} proposed different approaches for evaluating the three categories of losses and examined their effects on system responsiveness. $tan\varphi \u2032$, $tan\theta \u2032$, and $tan\delta \u2032$ are usually used to denote the intensive dissipation factors for elastic, piezoelectric, and dielectric materials, respectively. The parameter $\delta \u2032$ denotes the phase delay of the electric displacement to the applied electric field under the constant stress condition, $\varphi \u2032$ denotes the phase delay of the strain to the applied stress under the constant electric field condition, and $\theta \u2032$ is the phase delay of the strain to the applied stress under the constant electric field condition. In this study, the mechanical quality factor (denoted by $Qm$), which is the inverse of the intensive mechanical loss $(tan\varphi \u2032)$, is estimated for resonance and anti-resonance frequencies of distinct vibration modes using the half-power bandwidth approach.^{22} The piezoelectric loss is expected to have a negligible influence on the overall quality factor,^{5,23} and the dielectric quality factor $(Qd=1/tan\delta \u2032=750)$ was provided by the manufacturer. It has also been noticed that the mechanical loss dominates the in-air damping ratio.^{21}

Figure 2(a) depicts the resonance and anti-resonance frequencies of several piezoelectric disk vibration modes. Radial (R), edge (E), high-frequency radial (A), thickness-shear (TS), and thickness-extensional (TE) modes are among the symmetric modes (bilateral symmetry about the disk midline in the y-axis). The displacement profiles of the disk mode shapes are shown in Fig. 2(b) for better identification of the vibration modes associated with the TWUPT system. The letter and number represent the type and order of the mode, respectively. The disk operating at R modes features radial distortion patterns in the lateral direction and grows in the thickness direction due to Poisson's ratio effects,^{24} as shown in Fig. 2(b). The disk has a substantial axial displacement at the edge of the disk in the E mode. It is observed to be comparable to how the center of the disk mechanically clamps the edge. In particular, in the E mode, the disk's edge experienced significant axial movement, indicating that vibration energy is restricted to the disk's edge.

The shear displacement in the TS mode is parallel to the polling direction (thickness direction). The TS mode is defined by cylindrical rings inside the disk that are parallel to the thickness direction and have no radial displacement. For TE mode, the modes' axial displacements are similar to that of the radial modes as they both change with the radial position around a mean value. In the case of radial modes, this mean value is zero; hence, the nodal circles are the sites where the displacement has its mean value. The axial displacement pattern is identical to that of the radial modes but overlaid on a continuous shift in the case of the TE modes, where the mean value is non-zero. In conclusion, we compare the axial displacement of the piezoelectric disk surface of distinct vibration modes as a technique for identifying the modes. With respect to radial location, the axial displacements of vibration modes range from the most negative (compression) to the most positive (tension). The average axial displacement along the radial line of various modes is of interest for distinctively differentiating vibration modes. The radial modes' nodal circles represent the zero-mean value, but the thickness-extensional modes have a non-zero mean value. The high-frequency radial modes (A modes), as seen in Fig. 2(b), are the last set of modes. The axial surface displacement in the high-frequency modes (A) is essentially a superposition of the low-order R-mode displacement and the high-order R-mode displacement. On the other hand, the radial displacement at the cylindrical surface is found to be negligible. It should be noted that our investigation is focused on a frequency range of 0–600 kHz, which includes a variety of vibrational modes up to the first thickness-extensional mode, Fig. 2(a). Since thickness shear (TS) modes are absent in the study range, they are not shown in Fig. 2(b).

In Fig. 3, we plot the variations of the TWUPT system's electrical impedance with the excitation frequency. These data are used in a modal analysis process to better understand the dynamics associated with the TWUPT system and provide insight into its response limits. We note that FRFs were computed using a FEM after evaluating the mechanical and electromechanical material parameters of the transmitter and receiver, as well as the elastic wall and bonding layers, as reported in Table I. The FEM simulations and experimental measurements are carried out on the setup described in Fig. 1, which is connected to 100 Ω load. We observe that the response characteristics of the TWUPT system are fulfilled with additional vibration modes when the impedance of the PZT disk is satisfied, as shown in Fig. 3. The comparison of the experimental results with their FEM numerical counterparts, including the resonance and anti-resonance frequencies and amplitudes of the TWUPT system's various vibration modes under free-free boundary conditions, is also shown in Fig. 3. We can see a slight difference in the frequency and amplitude of the electrical impedance at the resonance frequencies of all modes in Figs. 2 and 3. This discrepancy can be explained in a variety of ways. It may have manifested as a result of the non-uniqueness of the material parameter modifications, as well as the lack of proportional damping.

In Fig. 3, each zoomed-in plot depicts the identified modes at the selected frequency range. Furthermore, the visualization of the mode shapes obtained from FEM simulations in Fig. 4 allows for more accurate mode identification. These modes include bilateral symmetric R, E, and TE as well as asymmetric flexural (F) modes. In the case of the single disk, the asymmetric modes (F modes) cannot be activated by the input voltage to the fully electroded parallel surfaces of the disk.^{25} Due to asymmetric mechanical boundary conditions, adding further layers to the piezoelectric disk, as shown in Fig. 1, results in a multiplex electromechanical structure with activated F vibrating modes, as illustrated in Figs. 3 and 4(a). F modes have radial symmetry about the z-axis but bilateral asymmetry about the y-axis at the disk midline. Local resonances associated with standing waves in the elastic layers, i.e., aluminum layer, shown in Fig. 4(b) are barely observed in the impedance plots in Fig. 3 due to the fact that, unlike piezoelectric elements, elastic layers are not electrically coupled.

### B. TWUPT model validation

The geometric and material characteristics provided in Fig. 1 and Table I are used to conduct numerical simulations using multiphysics FEM in COMSOL®. A free triangular mesh with an element size of $(1/10\lambda i)$, where $\lambda i$ represents the acoustic wavelength in the *i*th domain, is used to discretize the domains in a 2D axisymmetric acoustic-structure model. It is worthy to note that the proposed approach in this paper involves 3D optimization, with the tuned parameters being the radius, thickness, and rotation angle. The method is critical for asymmetric electrical excitation in piezoelectric transducer disks with wrap-around electrodes. Wrap-around electrodes have the electrode wrapped from the bottom to the top surface. In TWUPT systems, where electrical connections are formed on only one side, asymmetric excitation is achieved by using wrap-around-electrode piezoelectric disks. It is because of their ease of installation that the attached side can be bonded to the elastic wall. On the other hand, for the specific case of symmetric excitation, a 2D axisymmetric model is used to solve the physics of the TWUPT system at a lower computational cost. To ensure proper sampling of the wave at the selected range of frequencies employed in the study, a factor of ten is used. Experiments are carried out to excite the through-wall system with free-free boundary conditions (see Fig. 1) in order to validate the FEM simulations and get a better understanding of the coupled system dynamics. The KEYSIGHT 33500B function generator created continuous sinusoidal waves that are delivered to the transmitter at specific frequencies. The electrodes of the piezoelectric receiver disk are shunted to a resistance substitution box, IET labs RS-2W. The input voltage, $Vt$, input current (measured using Tektronix TCP2020), $It$, and output voltage, $Vr$, are recorded using a Tektronix TBS20000B digital oscilloscope with 1 GHz sampling frequency. The efficiency is defined as $\eta =Po/Pi$, where $Po=Vr2/Rl$ is the output power dissipated at the resistive load, $Rl$, attached to the receiver and $Pi=Re(12VtIt\u2217)$ is the input active power received by the transmitter. The asterisk denotes complex conjugate. *Re* denotes the real part of the complex number.

To determine the optimal load resistance, we use a gradient-based optimization method (Matlab built-in function fmincon) to maximize the system's efficiency and output power. This function enables solving for a local minimum of a constrained nonlinear multivariable function using a gradient-based technique. This optimization tool was selected in this study given its inexpensive processing costs and good accuracy when dealing with a small number of variables.^{5} When dealing with a large number of variables, the GA, which was used for parameter identification in Sec. II A, is more efficient. The FEM's objective functions and tolerance were all employed in conjunction with the optimization function. At an excitation frequency of 42.2 kHz (first radial mode “R1”) and an ideal load resistance of $Rl=456\Omega $, the highest efficiency of 79.9% is achieved. We test the systems' power efficiency and input-to-output voltage ratio FRFs for excitation frequencies spanning a wide range, including the first thickness extensional mode at 437.7 kHz at $470\Omega $ load resistance (experimental optimum load). We show in Fig. 5 a comparison between experimental and simulation results. A good agreement between the experimental data and their FEM numerical counterparts is obtained. The slight discrepancy is mostly associated with fabrication flaws producing nonuniform thickness of epoxy coupling layers.

As demonstrated in Figs. 5(a) and 5(c), the peaks in the normalized voltage and efficiency FRFs (TE1 and R modes) match with the TWUPT system's resonance frequencies shown in Fig. 3. Figure 5(b) is an enlarged image of Fig. 5(a) near the PZT disk's first radial mode at 40.5 kHz, which constitutes the most efficient working mode of the present system. The first radial (R1) and second flexural (F2) modes are clearly close to each other. In Sec. II C, the impact of such mode coupling on the overall system's performance will be examined.

Thickness mode devices have a low aspect ratio (diameter-to-thickness ratio, AR < 20) and relatively lower thickness mode resonance frequencies. For these transducers, the displacement fields associated with the TE mode become nonuniform, and the coupling between the disk's radial and high-frequency radial modes (A modes) with TE1 mode is more pronounced,^{25–27} as shown in Figs. 2 and 3 for AR = 12. This can have a negative impact on the transducer's operation. As a result, low aspect ratio transducers operating around the TE1 mode are not suitable for TWUPT systems. A pure radial mode free of interference can be achieved by carefully selecting the appropriate system dimensions. According to Fig. 5(c), the highest efficiency occurs in the vicinity of the first radial mode. This is because, unlike lower-order radial modes, the interference at the TE1 mode has a negative impact on the overall efficiency.

### C. Geometric-dependent mode couplings: Hills and valleys of efficiency in TWUPT systems

According to Fig. 5(c), the highest efficiency occurs in the vicinity of the first radial mode at 42.2 kHz. As a result, we consider a radial mode device and optimize the system's efficiency for various transmitter and receiver thicknesses. Figures 6, 7(a), and 7(b) show the FEM numerical results for fixed elastic wall thicknesses of $tw=6$ and $3mm$, respectively. Each point on the contour plot shown in the aforementioned figures represents an electrically optimized TWUPT system (optimized in terms of frequency and attached resistive load).

The contour plot shown in Fig. 6(a) displays the variations of the system's efficiency with transmitter and receiver thicknesses. It reveals a valley that extends the depression of the efficiency data, typically between hill ranges. This is explained by the occurrence of mode coupling. The first mode coupling occurs between the disks' first radial mode (R1, mode of operation) and the second flexural mode (F2) in a nearly quarter-ring pattern, as shown in Fig. 6(a). Figure 6(b) shows the voltage output FRFs normalized to the input voltage. The R1 and F2 modes are weakly coupled and merging in Fig. 6(b-1) which leads to an efficiency of 71.4% [point 1 in Fig. 6(a)]. The interaction between the two modes is complete in Fig. 6(b-2) which corresponds to a valley, resulting in an efficiency of 68.8% [point 2 in Fig. 6(a)]. Then, the two modes start to separate, as shown in Fig. 6(b-3), which corresponds to an efficiency of 77.7% [point 3 in Fig. 6(a)]. These two modes are finally well separated, as shown in Fig. 6(b-4), which corresponds to an efficiency of 79% [point 4 in Fig. 6(a)]. The R1 mode frequency of the TWUPT system with equal transmitter and receiver thicknesses $(tt=tr)$ calculated using FEM is shown in Fig. 6(c). There is a strong coupling between the first radial (R1) and second flexural (F2) modes for systems with roughly $3mm<(tt=tr)<4mm$, as can be seen, which corresponds to the valley observed in Fig. 6(a) denoted by (2). As seen in Fig. 6, there is a 1:1 mode coupling between the first symmetric R1 and second antisymmetric F2 modes through an internal resonance. It should be noticed that once the two modes intersect, the internal resonance properties predominate. The TWUPT system's R1 and F2 uncoupled vibration modes are also depicted in Fig. 4 for further understanding (a).

Figures 7(a) and 7(b) depict the effects of transmitter and receiver thicknesses on the system's efficiency for $tw=6and3mm$, respectively, for extended $tt$ and $tr$ values. The valley of efficiency data can be explained by an increase in the overall damping, which leads to increased energy loss.^{28} Another interference occurs between the first radial mode (R1, mode of operation), the first edge mode (E1), and the first flexural mode (F1), resulting in a stronger interference and thus lower efficiencies. In fact, the interference is caused primarily by the relative change in the frequencies of the flexural modes (F) to the first radial mode (R1) frequency when the system parameters are changed. Changing the wall thickness from $tw=6to3mm$ reduces the flexural frequencies, causing interference to occur in systems with thicker piezoelectric disks. This can be noticed in Fig. 7(b) as a larger radius ring pattern is obtained in comparison to the pattern shown in Fig. 7(a). The occurrence of mode coupling is also indicated by the abrupt change in the optimum load and operating frequencies, which are depicted as shadowed regions in Figs. 7(c) and 7(d).

It is worth noting that the losses are frequency-dependent. As a result, radial mode devices with thin piezoelectric elements have higher efficiencies than thickness mode devices. This is because, unlike thickness mode devices, varying the thickness of the piezoelectric disks has only a minor effect on the operating frequency (35–50 kHz), as shown in Figs. 7(c) and 7(d).

It is a well-known design criterion to drive the system's transducer at its resonance frequency for improved performance. The inclusion of losses in the TWUPT model, however, shifts the optimal operating frequency away from the resonance frequency. Operating at resonance generates a lot of heat, resulting in a low quality-factor, leading to a low efficiency. Figure 8 explains the abrupt decrease in operating frequency reported in Fig. 7(c). It depicts the impedance and efficiency curves of four different systems selected from the efficiency map in Fig. 7(a) as a function of frequency near the optimal operating frequency, which is marked by a dashed vertical line. A weak mode coupling occurs in Fig. 8(a), resulting in 74% efficiency. The optimal operating frequency shifted away from F2 mode (denoted by a black dashed line) and toward the anti-resonance frequency of R1 mode (green dashed line). Figure 8(b) depicts the onset of the region of low efficiencies where the R1 and F2 modes interact. At that point, the system impedance reaches its maximum and, thus, the matching load is also at its peak. Figure 8(c) depicts the point at which mode interference has the greatest impact on the system's efficiency. The operating point is at the radial resonant frequency. Finally, as depicted in Fig. 8(d), the efficiency recovers to 79.9% and the operating frequency is nearly midway between the resonance and antiresonance frequencies of the radial mode. Experimentally, this observation was confirmed in Fig. 8(d). As reported in Ref. 17, the maximum quality factor point, and thus the maximum efficiency, occurs somewhere between the resonance and antiresonance frequencies of a given vibration mode.

Next, we optimize the system's efficiency for various wall geometries. We show in the contour plot of Fig. 9 the variations of the optimized efficiency with the wall thickness and radius. We note that each point on the plot represents an electrically optimized system, as previously stated. Clearly, increasing the wall thickness results in a uniform decrease in efficiency with slight fluctuations. The increase in the total energy dissipation due to damping in the wall when increasing the wall thickness is mostly the cause of the observed reduction in the efficiency. On the other hand, the system's efficiency varies with the radius for a given wall thickness. Mode coupling with asymmetric modes can also lead to such behavior. Changing the wall radius causes the flexural mode frequencies to change at a faster rate than the radial mode frequencies, causing modes to approach and cross each other and multiple modes coupling to occur. In fact, changing the wall radius has only a slight effect on the radial mode (R1) frequency.

The results in Fig. 9 demonstrate that system 1 (same as depicted in Fig. 1 and analyzed in Figs. 2–8) is not the most efficient system due to the proximity of the second flexural mode (F2) to the first radial mode (R1), as previously demonstrated. Here, the results reveal an alternative method for avoiding mode coupling by changing the wall radius, thereby increasing the system's efficiency. Thus, to confirm the simulation results, we build a second TWUPT system (denoted by system 2 in Fig. 9) with the same dimensions but a wall radius of *r _{w}* = 42.5 mm. The efficiency of system 2 is found experimentally equal to 85%. This value agrees with the efficiency obtained from FEM simulations, which is found equal to 86.5%. Based on the results reported in Fig. 9, radial mode TWUPT devices with a larger radius wall than the radius of the piezoelectric transducer can exhibit weaker mode coupling with higher-order flexural modes and, thus, have a better chance to achieve higher efficiencies. The material damping, on the other hand, limits the ability to increase the wall radius.

In addition to the efficiency results shown in Fig. 7(a), we formulate a new optimization problem based on maximizing the output active power. The corresponding results are depicted in Fig. 10. It is important to note that the systems analyzed in Figs. 7(a) and 10 operate at different frequency and resistive load values, which means that unless re-tuned, the specific system in Fig. 10 will not have the efficiency shown in Fig. 7(a). For different receiver thicknesses up to about 3 mm, the result predicts a high output power per unit square of input voltage when the transmitter has a very small thickness. In order to achieve high output power levels for a given input voltage, two main conditions must be met: high input active power and a high-efficiency system. Higher input active power values can be obtained with higher input current and power factor levels. Higher current values can be obtained with higher piezoelectric capacitance values found in thin piezoelectric elements. Power factor is defined as the ratio of real or active (dissipated) power to apparent (available) power and it is a measure of how in phase the input current is with the voltage. Figure 10 mainly shows how the maximum capacitance and power factor points approach each other where the highest power values appear (the effects of transmitter and receiver thicknesses on the input power factor were reported in Fig. S2 of the supplementary material). It is also observed that mode coupling occurrence has no clear effect on the optimized output power. It is worth mentioning that an external series or parallel inductor is commonly used to compensate for the TWUPT system's capacitance, hence increasing the power factor leading to higher output power levels.^{29}

### D. Dependence of the electromechanical coupling factor on the vibration modes

The electromechanical coupling factor is the ratio of electrical and mechanical energy exchanged throughout a piezoelectric element's excitation cycle. Methods for achieving a standard definition of the electromechanical coupling factor have been devised in several studies. Some of these techniques are centered on energy efficiency.^{30} Under dynamic conditions, the dynamic electromechanical coupling factor, $kd$, of a piezoelectric element associated with the internal energies is defined in Eq. (4).^{30} The elastic, dielectric, and mutual piezoelectric energies are represented by the piezoelectric internal energy components $U1$, $U2$, and $U12$, respectively, in Eqs. (1)–(3). It is clearly noted that the associated energy densities are integrated over the volume of the piezoelectric disk in order to account for the distribution of stress and electric field. Because of the phase differences between the stress and electric field owing to material damping, the resulting energy is a complex number. As stated in Eq. (4),^{30} the variable $kd$ is defined by considering the real part of the computed energy components in Eqs. (1)–(3),

where $sijE$, $\epsilon ijT$, and $dij$ denote the elastic compliance under constant electric field *E*, the dielectric constant under constant stress *T*, and the piezoelectric constant, respectively. The constants are obtained from Table I. The variables $Ti$ and $Ei$ are the components of stress and electric field vectors, respectively. The asterisk denotes complex conjugate and *Re* denotes the real part of the complex number.

As shown in Figs. 11(a)–11(d), we report the variation of the dynamic coupling factor $kd$ with excitation frequency and piezoelectric layer thicknesses around the R1 mode (mode of operation) under open-circuit conditions (R → ∞). The arrows added to the contour plots show the change in resonance frequencies of the R1 and F2 modes as the thicknesses of the piezoelectric layers are increased. Both modes are clearly visible when the two arrows cross each other at $tt=tr\u22484.25mm$ and $tt=tr\u22485.5mm$ for the systems in Figs. 11(a) and 11(b), respectively. This interaction causes a slight decrease in the $kd$ value, as shown in Figs. 11(c) and 11(d). The exchange of energy between the modes as they interact can explain such behavior.

## III. EFFICIENT TWUPT SYSTEM BASED ON MODE COUPLINGS

### A. Trade-off between efficiency and output power

Previous results demonstrated that the system's operating points (frequency and attached load) to achieve maximum efficiency and output power differ. In fact, maximum efficiency and power cannot be obtained simultaneously. For system 1, a maximum efficiency of 79.9% corresponds to an output power of $0.88mW/Vt2$ and a maximum power of $2.2mW/Vt2$ corresponds to an efficiency of 49%. These correspond to the two extreme points highlighted in Fig. 12 by circles. Such a trade-off between efficiency and power is investigated. Figure 12 depicts the maximum output power of two different TWUPT systems for a given efficiency. We observe that the efficiency is inversely proportional to the output power. System 2 constitutes the optimized system because it avoids mode coupling, as previously mentioned, and thus produces higher efficiency levels than system 1 for each given output power unit. Depending on the application, the system can be tuned to achieve a medium efficiency and power output. By setting the excitation frequency equal to 42.2 kHz and attaching a load of 1570 Ω, for example, a moderate efficiency of 70% and an output power of $1.74mW/Vt2$ can be obtained.

### B. Effects of circuit load and frequency on the overall efficiency

The optimum system (system 2) has a matching load of 650 Ω and an operating frequency near R1 at 40 kHz. According to the numerical simulations, tuning the system at the given frequency and load results in an efficiency of 86.5%. The experiments, in agreement, reveal a maximum efficiency of around 85% at a frequency of 41 kHz and a resistive load of 750 Ω. The experimental findings corroborate the simulations. In the optimal configuration, there is a slight frequency shift. The efficiency varies with the excitation frequency and resistive load, as seen in Fig. 13. The results are obtained from the experiments and the FEM simulations. Similar features are obtained. This demonstrates the capability of the FEM to properly capture the inherent physical aspects of the system. The numerical and experimental results reveal that varying the resistive load in the vicinity of the optimal value has an insignificant effect on the system's efficiency as shown in Fig. 13. As a result, a high tolerance for selecting the attached resistive load is not required.

## IV. CONCLUSIONS

The effects of elastic wave propagation on mode couplings in multiplex electromechanical structures, such as multi-layered TWUPT systems, were studied. We investigated mode couplings in an optimal scenario, including their physical origins, modeling, and weak and strong coupling regimes. The system layout's size (volume), operating frequency, and matching circuit load optimization were all taken into account. The use of a GA for system characterization, as well as a gradient-based optimization tool in conjunction with FEM simulations to maximize efficiency and electrical output power, were highlighted and verified with the experiments. To demonstrate how eigenmode features affect optimization results, several numerical and experimental studies were done. The interaction of symmetric radial and asymmetric flexural modes has been found to enhance damping and consequently reduce the overall efficiency of the device.

We developed a three-dimensional FEM that includes system losses for optimizing TWUPT systems. Our model enables us to assess the feasibility of efficiently running systems with a wide range of transducer aspect ratios in various vibration modes. However, the previously proposed one-dimensional modeling methodologies, such as analytical and equivalent circuit modeling, could not be used to conduct such a study. According to numerical calculations and experiments, the optimal system was found at a frequency approaching R1 mode, with an efficiency of 86.5%. Under open-circuit conditions, the dynamic electromechanical coupling factor was shown to vary with excitation frequency and piezoelectric layer thicknesses around the R1 mode (mode of operation). As the thicknesses of the piezoelectric layers are increased, the resonance frequencies of the R1 and F2 modes shift. The interplay between the modes leads the electromechanical coupling factor value to drop slightly. Such a phenomenon can be explained by the flow of energy between the modes as they interact.

Based on the acquired results, the primary design factors that enable to maximize overall efficiency can be summarized: (1) Use of radial mode devices to achieve high-efficiency levels. In fact, as the aspect ratio of the transducer decreases due to the inevitable mode coupling near the thickness mode, the efficiency gap between thickness mode and radial mode devices widens. The piezoelectric thicknesses need to be selected so that no mode coupling occurs with the radial mode. The operation near the radial mode leads to superior performance because the first radial mode occurs at a lower frequency than that the thickness mode and then the system experiences fewer frequency-dependent losses. Furthermore, unlike the thickness mode, asymmetric mode coupling with the radial mode can be easily avoided for all aspect ratios. This type of coupling is unavoidable, but it can be overlooked in thickness-mode devices with high-aspect-ratio transducers (typically AR > 20). (2) Avoid mode coupling of asymmetric modes (F modes) to the mode of operation because it increases the damping and, thus, decreases the overall efficiency. (3) Avoid radial mode systems with very thin receivers. (4) Use equal-diameter transmitter and receiver radial mode systems. This is because the resonance frequencies of the transmitter and receiver are matched.

Because efficiency and power optimization cannot be achieved concurrently for a given system under excitation, another power optimization study was conducted. A trade-off analysis was also carried out between output power and efficiency. It demonstrated the possibility of sacrificing some efficiency in exchange for moderate output power, which can be useful for a variety of applications such as wireless ultrasound power and data delivery systems.

## SUPPLEMENTARY MATERIAL

See the supplementary material for a more detailed analysis and numerical optimization results representing frequency-varying material properties that were used in this study. The effects of the thickness of the transmitter and receiver on the input power factor are also demonstrated.

**ACKNOWLEDGMENTS**

M.S.A. and S.S. gratefully acknowledge support from U.S. National Science Foundation (NSF) under Grant No. ECCS-1711139 and Institute for Critical Technology and Applied Science (ICTAS) at Virginia Tech.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Moustafa Sayed Ahmed:** Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **Mehdi Ghommem:** Conceptualization (equal); Supervision (equal); Writing – review & editing (equal). **Shima Shahab:** Conceptualization (equal); Funding acquisition (lead); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: PIEZOELECTRIC CONSTANTS AND LINEAR CONSTITUTIVE EQUATIONS

## REFERENCES

_{2}nanoelectromechanical system