# Mapped orthogonal functions method applied to acoustic waves-based devices

Jean-Etienne Lefebvre, J.G. Yu, Faniry Emilson Ratolojanahary, Lahoucine

Elmaimouni, Wei-Jiang Xu, Tadeusz Gryba

► To cite this version:

Jean-Etienne Lefebvre, J.G. Yu, Faniry Emilson Ratolojanahary, Lahoucine Elmaimouni, Wei-Jiang Xu, et al.. Mapped orthogonal functions method applied to acoustic waves-based devices. AIP Advances, 2016, 6 (6), 065307, 18 p. 10.1063/1.4953847. [hal-04080549](https://hal-04080549)

HAL Id: hal-04080549

https://uphf.hal.science/hal-04080549

Submitted on 13 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Mapped orthogonal functions method applied to acoustic waves-based devices

Jean-Etienne Lefebvre, J.G. Yu, Faniry Emilson Ratolojanahary, Lahoucine

Elmaimouni, Wei-Jiang Xu, Tadeusz Gryba

► To cite this version:

Jean-Etienne Lefebvre, J.G. Yu, Faniry Emilson Ratolojanahary, Lahoucine Elmaimouni, Wei-Jiang Xu, et al.. Mapped orthogonal functions method applied to acoustic waves-based devices. AIP Advances, 2016, 6 (6), 065307, 18 p. 10.1063/1.4953847. hal-04080549

HAL Id: hal-04080549

https://uphf.hal.science/hal-04080549v1

Submitted on 13 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mapped orthogonal functions method applied to acoustic waves-based devices

J. E. Lefebvre 🐦; J. G. Yu; F. E. Ratolojanahary 🐦; L. Elmaimouni; W. J. Xu; T. Gryba

<figure><img src="image_2.png" /></figure>

Check for updates

AIP Advances 6, 065307 (2016)

https://doi.org/10.1063/1.4953847

<figure><img src="image_3.png" /></figure>

View
Online

<figure><img src="image_4.png" /></figure>

Export
Citation

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>
# Mapped orthogonal functions method applied to acoustic waves-based devices

J. E. Lefebvre 🐦; J. G. Yu; F. E. Ratolojanahary 🐦; L. Elmaimouni; W. J. Xu; T. Gryba

<figure><img src="image_3.png" /></figure>

Check for updates

AIP Advances 6, 065307 (2016)

https://doi.org/10.1063/1.4953847

<figure><img src="image_4.png" /></figure>

View
Online

<figure><img src="image_5.png" /></figure>

Export
Citation

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>
# Mapped orthogonal functions method applied to acoustic waves-based devices

J. E. Lefebvre,¹,ª J. G. Yu,² F. E. Ratolojanahary,³ L. Elmaimouni,⁴ W. J. Xu,¹ and T. Gryba¹

¹UVHC, IEMN-DOAE, 59313 Valenciennes Cedex 9, France; CNRS, UMR 8520, F-59650 Villeneuve d'Ascq, France

²College of Aerospace Engineering, Chongqing University, P.R. China, 400044

³LAPAUF, Université de Fianarantsoa, 301 Fianarantsoa, Madagascar

⁴LSIE-ERMAM, Faculté Polydisciplinaire d'Ouarzazate, Université Ibn Zohr, 45000 Ouarzazate, Morocco

(Received 26 February 2016; accepted 30 May 2016; published online 8 June 2016)

This work presents the modelling of acoustic wave-based devices of various geometries through a mapped orthogonal functions method. A specificity of the method, namely the automatic incorporation of boundary conditions into equations of motion through position-dependent physical constants, is presented in detail. Formulations are given for two classes of problems: (i) problems with guided mode propagation and (ii) problems with stationary waves. The method's interest is demonstrated by several examples, a seven-layered plate, a 2D rectangular resonator and a 3D cylindrical resonator, showing how it is easy to obtain either dispersion curves and field profiles for devices with guided mode propagation or electrical response for devices with stationary waves. Extensions and possible further developments are also given. © 2016 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). [http://dx.doi.org/10.1063/1.4953847]

## I. INTRODUCTION

Acoustic waves play a key role in many fields such as non-destructive testing (NDT) for material inspection or evaluation and also, with the advances in reliability, structural health monitoring (SHM) with continuous screening and early warning of defect growth and structural failure. Some applications use bulk waves which cover only a small localized section of a structure but a majority of applications rely on guided waves which provide a complete coverage of the waveguide cross-section. The applications are in aircraft, composite materials, pipeline and tank inspection.

All those applications encompass a great variety of geometries, plates, rods, hollow cylinders, spheres, multilayer structures, layer or multiple layers on a half-space and curved surfaces. Accurate and efficient numerical methods are needed to calculate the fundamentals of wave propagation in such geometries, guided wave dispersion curves and associated field profiles. Among computational methods that have been used are the transfer matrix approach developed by Thomson¹ and the global matrix method by Knopoff² applicable to multilayered plates and cylindrical structures.³ The global matrix method was proposed to overcome the problem of numerical instability of the solution given by the transfer matrix method, known as the large *fd* problem when layers of large thickness *d* are present and high frequencies *f* are being considered. This is crucial when inhomogeneous waves are considered. The global matrix method is robust, it does not suffer from numerical instability, but is computationally expensive in time and memory resources.

Later on, the transfer matrix method was improved to account for anisotropic materials and damping.⁴,⁵ Also a new alternative method to overcome the numerical instability problem was

<footnote>ªjean-etienne.lefebvre@univ-valenciennes.fr</footnote>
# Mapped orthogonal functions method applied to acoustic waves-based devices

J. E. Lefebvre,¹,ª J. G. Yu,² F. E. Ratolojanahary,³ L. Elmaimouni,⁴ W. J. Xu,¹ and T. Gryba¹

¹UVHC, IEMN-DOAE, 59313 Valenciennes Cedex 9, France; CNRS, UMR 8520, F-59650 Villeneuve d'Ascq, France

²College of Aerospace Engineering, Chongqing University, P.R. China, 400044

³LAPAUF, Université de Fianarantsoa, 301 Fianarantsoa, Madagascar

⁴LSIE-ERMAM, Faculté Polydisciplinaire d'Ouarzazate, Université Ibn Zohr, 45000 Ouarzazate, Morocco

(Received 26 February 2016; accepted 30 May 2016; published online 8 June 2016)

This work presents the modelling of acoustic wave-based devices of various geometries through a mapped orthogonal functions method. A specificity of the method, namely the automatic incorporation of boundary conditions into equations of motion through position-dependent physical constants, is presented in detail. Formulations are given for two classes of problems: (i) problems with guided mode propagation and (ii) problems with stationary waves. The method's interest is demonstrated by several examples, a seven-layered plate, a 2D rectangular resonator and a 3D cylindrical resonator, showing how it is easy to obtain either dispersion curves and field profiles for devices with guided mode propagation or electrical response for devices with stationary waves. Extensions and possible further developments are also given. © 2016 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). [http://dx.doi.org/10.1063/1.4953847]

## I. INTRODUCTION

Acoustic waves play a key role in many fields such as non-destructive testing (NDT) for material inspection or evaluation and also, with the advances in reliability, structural health monitoring (SHM) with continuous screening and early warning of defect growth and structural failure. Some applications use bulk waves which cover only a small localized section of a structure but a majority of applications rely on guided waves which provide a complete coverage of the waveguide cross-section. The applications are in aircraft, composite materials, pipeline and tank inspection.

All those applications encompass a great variety of geometries, plates, rods, hollow cylinders, spheres, multilayer structures, layer or multiple layers on a half-space and curved surfaces. Accurate and efficient numerical methods are needed to calculate the fundamentals of wave propagation in such geometries, guided wave dispersion curves and associated field profiles. Among computational methods that have been used are the transfer matrix approach developed by Thomson¹ and the global matrix method by Knopoff² applicable to multilayered plates and cylindrical structures.³ The global matrix method was proposed to overcome the problem of numerical instability of the solution given by the transfer matrix method, known as the large *fd* problem when layers of large thickness *d* are present and high frequencies *f* are being considered. This is crucial when inhomogeneous waves are considered. The global matrix method is robust, it does not suffer from numerical instability, but is computationally expensive in time and memory resources.

Later on, the transfer matrix method was improved to account for anisotropic materials and damping.⁴,⁵ Also a new alternative method to overcome the numerical instability problem was

<footnote>ªjean-etienne.lefebvre@univ-valenciennes.fr</footnote>
proposed which consists in the introduction of a layer stiffness matrix associated with an efficient
recursive algorithm to calculate the global stiffness matrix for an arbitrary anisotropic layered
structure.⁶ A recursive surface impedance method was also proposed.⁷

With the matrix methods, a solution requires the determination of the roots of a characteristic
function established from equations of motion and continuity and boundary conditions. Efficient
numerical iteration techniques have to be employed to perform this search.³,⁸ With leaky waves or
with damped media, the wavenumber and the characteristic function become complex and the roots
search is very challenging. Alternatives to the classical solution of a characteristic function have
been proposed.⁹⁻¹¹ They allow to compute real and complex wavenumbers via a classical matricial
eigenvalue problem obtained from the motion equations and the continuity and boundary conditions
along with the use of a polynomial expansion approach.

Functionally graded materials (FGM) are a specific class of heterogeneous materials. With continuously varying physical properties, these materials exhibit interesting mechanical and/or thermal behavior and are developed in nature (wood, bone...) and industry processes. Transversely continuously inhomogeneous waveguides can be dealt with thanks to a matrix method by approximating the waveguide material by successive layers of homogeneous materials whose properties vary slowly from layer to layer and whose thickness is small compared to wavelength. A major drawback of the process is that it involves, by modelling the FGM by homogeneous layers, an approximation before any calculation. An alternative method consists in using geometrical integrators in the form of a Magnus expansion¹² to compute the transfer matrix for a thin vertically nonhomogeneous layer and an efficient recursive algorithm to compute from the thin layer solution the wave propagation solution for an arbitrary FGM layered structure. For an overview on approximate methods appropriate to FGM, see Ref. 12 and references herein. Two alternative approaches without using a multilayered model for the FGM were proposed. The first one takes advantage of the knowledge of an analytical solution of the wave equation, the matricant, which is explicitly expressed via the Peano series expansion¹³⁻¹⁵ and uses a general method based on the Stroh's formalism.¹⁶ A specific advantage of the method is that it allows to evaluate the degree of round-off error and the validity domain. The second one uses polynomials to express spatially varying properties and elastic fields and a mapped orthogonal functions technique.¹⁷⁻²⁰ The latter was also extended to waveguides of arbitrary cross-sections.²¹

Alternative methods based on finite elements have also been developed to compute the fundaments of wave propagation in heterogeneous and inhomogeneous waveguides of arbitrary cross-section. The most popular is the semi-analytical finite element (SAFE) method which requires the finite element discretization of the cross-section of the waveguide and assumes harmonic motion along the wave propagation direction.²²,²³ An interesting review on SAFE method can be found in Ref. 24 where various improvements given to the initial SAFE method to deal with specific cases are addressed, for instance non propagative modes, leaky modes and damped modes.

So far, we only considered ultrasonic devices which can be regarded as of unbounded extent in
at least one spatial dimension and involve propagative waves. But acoustic waves also play a key
role for ultrasonic devices of finite extent in one or several spatial dimensions and which involve sta-
tionary waves. Applications are numerous with the development of microelectromechanical systems
(MEMS) which allow miniaturized elements and improved performance. We can cite resonators,
oscillators, transformers, filters, biological sensors, devices for control and manipulation of particles
in suspensions or droplets in emulsions.

Surface acoustic waves (SAW) and bulk acoustic waves (BAW) are widely used for all these
applications. Telecommunication industry with the tremendous growth of mobile communications
makes extensive use of resonators and filters with GHz operating frequency range. Numerous
models have been developed to predict the behaviour of SAW and BAW devices in response to
external loads: P-matrix model<sup>25</sup> and coupling-of-modes<sup>26</sup> (COM) for SAW devices, matrix model<sup>27</sup>
for BAW devices and for both of them equivalent circuit model,<sup>28,29</sup> mapped orthogonal functions
model<sup>30,31</sup> and finite element method (FEM).<sup>32,33</sup> The FEM is always interesting because it can
handle complicated geometries and electrode patterns of real device structures but can sometimes be
computationally intensive.
proposed which consists in the introduction of a layer stiffness matrix associated with an efficient recursive algorithm to calculate the global stiffness matrix for an arbitrary anisotropic layered structure.⁶ A recursive surface impedance method was also proposed.⁷

With the matrix methods, a solution requires the determination of the roots of a characteristic function established from equations of motion and continuity and boundary conditions. Efficient numerical iteration techniques have to be employed to perform this search.³,⁸ With leaky waves or with damped media, the wavenumber and the characteristic function become complex and the roots search is very challenging. Alternatives to the classical solution of a characteristic function have been proposed.⁹⁻¹¹ They allow to compute real and complex wavenumbers via a classical matricial eigenvalue problem obtained from the motion equations and the continuity and boundary conditions along with the use of a polynomial expansion approach.

Functionally graded materials (FGM) are a specific class of heterogeneous materials. With continuously varying physical properties, these materials exhibit interesting mechanical and/or thermal behavior and are developed in nature (wood, bone...) and industry processes. Transversely continuously inhomogeneous waveguides can be dealt with thanks to a matrix method by approximating the waveguide material by successive layers of homogeneous materials whose properties vary slowly from layer to layer and whose thickness is small compared to wavelength. A major drawback of the process is that it involves, by modelling the FGM by homogeneous layers, an approximation before any calculation. An alternative method consists in using geometrical integrators in the form of a Magnus expansion¹² to compute the transfer matrix for a thin vertically nonhomogeneous layer and an efficient recursive algorithm to compute from the thin layer solution the wave propagation solution for an arbitrary FGM layered structure. For an overview on approximate methods appropriate to FGM, see Ref. 12 and references herein. Two alternative approaches without using a multilayered model for the FGM were proposed. The first one takes advantage of the knowledge of an analytical solution of the wave equation, the matricant, which is explicitly expressed via the Peano series expansion¹³⁻¹⁵ and uses a general method based on the Stroh's formalism.¹⁶ A specific advantage of the method is that it allows to evaluate the degree of round-off error and the validity domain. The second one uses polynomials to express spatially varying properties and elastic fields and a mapped orthogonal functions technique.¹⁷⁻²⁰ The latter was also extended to waveguides of arbitrary cross-sections.²¹

Alternative methods based on finite elements have also been developed to compute the fundamentals of wave propagation in heterogeneous and inhomogeneous waveguides of arbitrary cross-section. The most popular is the semi-analytical finite element (SAFE) method which requires the finite element discretization of the cross-section of the waveguide and assumes harmonic motion along the wave propagation direction.²²,²³ An interesting review on SAFE method can be found in Ref. 24 where various improvements given to the initial SAFE method to deal with specific cases are addressed, for instance non propagative modes, leaky modes and damped modes.

So far, we only considered ultrasonic devices which can be regarded as of unbounded extent in at least one spatial dimension and involve propagative waves. But acoustic waves also play a key role for ultrasonic devices of finite extent in one or several spatial dimensions and which involve stationary waves. Applications are numerous with the development of microelectromechanical systems (MEMS) which allow miniaturized elements and improved performance. We can cite resonators, oscillators, transformers, filters, biological sensors, devices for control and manipulation of particles in suspensions or droplets in emulsions.

Surface acoustic waves (SAW) and bulk acoustic waves (BAW) are widely used for all these applications. Telecommunication industry with the tremendous growth of mobile communications makes extensive use of resonators and filters with GHz operating frequency range. Numerous models have been developed to predict the behaviour of SAW and BAW devices in response to external loads: P-matrix model²⁵ and coupling-of-modes²⁶ (COM) for SAW devices, matrix model²⁷ for BAW devices and for both of them equivalent circuit model,²⁸,²⁹ mapped orthogonal functions model³⁰,³¹ and finite element method (FEM).³²,³³ The FEM is always interesting because it can handle complicated geometries and electrode patterns of real device structures but can sometimes be computationally intensive.
This paper aims at explaining the fundamentals of the mapped orthogonal functions method and reviewing all its potentialities in characterizing the acoustic wave-based devices with guided mode propagation and with stationary waves.

In section II, it is shown how, for the mapped orthogonal functions method, boundary conditions are automatically incorporated into equations of motion through the device of position-dependent physical constants. In section III, the method of solution is described for acoustic wave-based devices first with guided mode propagation and second with stationary waves. In section IV, the method is applied to obtain first the dispersion curves in a seven-layered plate and second resonance and antiresonance frequencies along with electrical impedances of rectangular and cylindrical resonators. In section V, we discuss some extensions and address possible further developments of the method.

## II. BOUNDARY CONDITIONS

A specificity of the mapped orthogonal functions method is to automatically incorporate the mechanical and electrical boundary conditions into the equations of motion through position-dependent physical constants rather than having recourse to a characteristic equation. The theoretical bases of this specificity have never been presented. Datta and Hunsinger³⁴ just mentioned that the position-dependent physical constants, when substituted in field equations, lead through derivatives to delta functions multiplying the normal-stress components, thus ensuring that the normal stress is zero at the boundary.

The theoretical bases are detailed here for mechanical and electrical boundary conditions and it is shown how position-dependent physical constants can fulfil this role.

### A. Stress-free boundary condition

The translational equation of motion for a vibrating medium is given by

$$ \nabla \cdot \mathbf{T}(\mathbf{r}, t) = \rho(\mathbf{r}) \frac{\partial^2 \mathbf{u}(\mathbf{r}, t)}{\partial t^2}, \qquad (1) $$

where $\mathbf{T}$, $\rho$ and $\mathbf{u}$ denote respectively stress, mass density and mechanical displacement evaluated at the time $t$ at the point $P$ given by the position vector $\mathbf{r} = (r_1, r_2, r_3)$. The left-hand side of Eq. (1) represents the divergence of stress defined by the relation³⁵

$$ \nabla \cdot \mathbf{T}(\mathbf{r}, t) = \lim_{\delta V \to 0} \frac{\oint_{\delta V} \mathbf{T}(\mathbf{r}, t) \cdot \mathbf{n} \, dS}{\delta V}, \qquad (2) $$

where $\delta V$ is an elemental volume in some orthogonal coordinate system and $\delta S$ is the surface enclosing this volume $\delta V$. $\mathbf{n}$ denotes an outward-directed unit vector.

If we study acoustic vibrations in a material body occupying a volume $V$ with a stress-free bounding surface $S$, Eq. (1) can, for the purpose of incorporating the stress-free mechanical boundary condition into the equations of motion, be written

$$ \nabla \cdot [f(\mathbf{r}) \mathbf{T}(\mathbf{r}, t)] = \rho(\mathbf{r}) f(\mathbf{r}) \frac{\partial^2 \mathbf{u}(\mathbf{r}, t)}{\partial t^2}, \qquad (3) $$

where

$$ f(\mathbf{r}) = \begin{cases} 1 & \text{if } P \text{ is within the body, bounding surface included} \\ 0 & \text{if } P \text{ is outside the body} \end{cases}. $$

The function $f(\mathbf{r})$ describes the region occupied by the material body. Given a well known vector identity, Eq. (3) can be written

$$ [\nabla f(\mathbf{r})] \cdot \mathbf{T}(\mathbf{r}, t) + f(\mathbf{r}) [\nabla \cdot \mathbf{T}(\mathbf{r}, t)] = \rho(\mathbf{r}) f(\mathbf{r}) \frac{\partial^2 \mathbf{u}(\mathbf{r}, t)}{\partial t^2}, \qquad (4) $$
This paper aims at explaining the fundamentals of the mapped orthogonal functions method and reviewing all its potentialities in characterizing the acoustic wave-based devices with guided mode propagation and with stationary waves.

In section II, it is shown how, for the mapped orthogonal functions method, boundary conditions are automatically incorporated into equations of motion through the device of position-dependent physical constants. In section III, the method of solution is described for acoustic wave-based devices first with guided mode propagation and second with stationary waves. In section IV, the method is applied to obtain first the dispersion curves in a seven-layered plate and second resonance and antiresonance frequencies along with electrical impedances of rectangular and cylindrical resonators. In section V, we discuss some extensions and address possible further developments of the method.

## II. BOUNDARY CONDITIONS

A specificity of the mapped orthogonal functions method is to automatically incorporate the mechanical and electrical boundary conditions into the equations of motion through position-dependent physical constants rather than having recourse to a characteristic equation. The theoretical bases of this specificity have never been presented. Datta and Hunsinger³⁴ just mentioned that the position-dependent physical constants, when substituted in field equations, lead through derivatives to delta functions multiplying the normal-stress components, thus ensuring that the normal stress is zero at the boundary.

The theoretical bases are detailed here for mechanical and electrical boundary conditions and it is shown how position-dependent physical constants can fulfil this role.

### A. Stress-free boundary condition

The translational equation of motion for a vibrating medium is given by

$$ \nabla \cdot \mathbf{T}(\mathbf{r},t) = \rho(\mathbf{r}) \frac{\partial^2 \mathbf{u}(\mathbf{r},t)}{\partial t^2}, \quad (1) $$

where $\mathbf{T}$, $\rho$ and $\mathbf{u}$ denote respectively stress, mass density and mechanical displacement evaluated at the time $t$ at the point $P$ given by the position vector $\mathbf{r} = (r_1, r_2, r_3)$. The left-hand side of Eq. (1) represents the divergence of stress defined by the relation³⁵

$$ \nabla \cdot \mathbf{T}(\mathbf{r},t) = \lim_{\delta V \to 0} \frac{\oint_{\delta V} \mathbf{T}(\mathbf{r},t) \cdot \mathbf{n} \, dS}{\delta V}, \quad (2) $$

where $\delta V$ is an elemental volume in some orthogonal coordinate system and $\delta S$ is the surface enclosing this volume $\delta V$. $\mathbf{n}$ denotes an outward-directed unit vector.

If we study acoustic vibrations in a material body occupying a volume $V$ with a stress-free bounding surface $S$, Eq. (1) can, for the purpose of incorporating the stress-free mechanical boundary condition into the equations of motion, be written

$$ \nabla \cdot [f(\mathbf{r}) \mathbf{T}(\mathbf{r},t)] = \rho(\mathbf{r}) f(\mathbf{r}) \frac{\partial^2 \mathbf{u}(\mathbf{r},t)}{\partial t^2}, \quad (3) $$

where

$$ f(\mathbf{r}) = \begin{cases} 1 & \text{if } P \text{ is within the body, bounding surface included} \\ 0 & \text{if } P \text{ is outside the body} \end{cases}. $$

The function $f(\mathbf{r})$ describes the region occupied by the material body. Given a well known vector identity, Eq. (3) can be written

$$ [\nabla f(\mathbf{r})] \cdot \mathbf{T}(\mathbf{r},t) + f(\mathbf{r}) [\nabla \cdot \mathbf{T}(\mathbf{r},t)] = \rho(\mathbf{r}) f(\mathbf{r}) \frac{\partial^2 \mathbf{u}(\mathbf{r},t)}{\partial t^2}, \quad (4) $$
where with regard to the definition of $f(\mathbf{r})$

$$ \nabla f(\mathbf{r}) = -\mathbf{n} \delta_s, \qquad (5) $$

with $\delta_s$ the Dirac delta-function.

For any point inside the volume $V$, bounding surface excluded, Eq. (4) restitutes with $f(\mathbf{r}) = 1$ and $\nabla f(\mathbf{r}) = 0$ the usual equation of motion given in Eq. (1). To see how the $f(\mathbf{r})$ function ensures the stress-free boundary condition on the outer surface, it is necessary to come back to the divergence of stress. Let us evaluate it on an elementary box of volume $V_{\text{box}}$ as shown in Fig. 1.

The top surface of the box $S_t$ belongs to the stress-free surface and the box thickness is infinitesimal. $S_s$ and $S_b$ denote respectively the side and bottom surfaces. Let us evaluate Eq. (4) with integration over the volume of the box

$$ \iiint_{V_{\text{box}}} \{[\nabla f(\mathbf{r})] \cdot \mathbf{T}(\mathbf{r},t) + f(\mathbf{r}) [\nabla \cdot \mathbf{T}(\mathbf{r},t)]\} dV = \iiint_{V_{\text{box}}} \rho(\mathbf{r}) f(\mathbf{r}) \frac{\partial^2 \mathbf{u}(\mathbf{r},t)}{\partial t^2} dV. \quad (6) $$

The first term of the left-hand side, by taking advantage of $\delta_s$, gives

$$ \iiint_{V_{\text{box}}} [\nabla f(\mathbf{r})] \cdot \mathbf{T}(\mathbf{r},t) dV = \iiint_{V_{\text{box}}} -\mathbf{n} \delta_s \cdot \mathbf{T}(\mathbf{r},t) dV = \iint_{S_t} -\mathbf{n} \mathbf{T}(\mathbf{r},t) dS $$

and the second term by applying the divergence theorem

$$ \iiint_{V_{\text{box}}} f(\mathbf{r}) [\nabla \cdot \mathbf{T}(\mathbf{r},t)] dV = \iint_{S_t} \mathbf{n} \mathbf{T}(\mathbf{r},t) dS + \iint_{S_b} -\mathbf{n} \mathbf{T}(\mathbf{r},t) dS + \iint_{S_s} \mathbf{T}(\mathbf{r},t) dS. $$

The contribution from the side of the box is negligible since its height is infinitesimal. The right-hand side of Eq. (6) would be zero since the volume of the box is infinitely small. Thus Eq. (6) leads to

$$ \iint_{S_b} -\mathbf{n} \mathbf{T}(\mathbf{r},t) dS = 0. \qquad (7) $$

Assuming that the cross-sectional area of the box is small enough that the stress is constant across the bottom surface then the normal stress vanishes on and near the bounding surface.

In conclusion, thanks to the function $f(\mathbf{r})$, Eq. (3) incorporates both Newton's second law Eq. (1) and the stress-free boundary condition.

## B. Electrically free boundary condition

For piezoelectric materials, the same device can be applied with the Maxwell equation

$$ \nabla \cdot \mathbf{D} = 0 \qquad (8) $$

to impose, when the electric fields outside the body are assumed to vanish, the electric boundary condition $D_n = 0$ for the normal electric displacement $D_n$ at the electrically free bounding surface.

In theory the electric fields are nonzero outside the body but the error introduced by this assumption is small because most piezoelectric materials have a high dielectric constant.²² For

<figure><img src="image_8.png" /><figcaption>FIG. 1. Elementary box for the evaluation of the divergence of stress.</figcaption></figure>
where with regard to the definition of $f(\mathbf{r})$

$$ \nabla f(\mathbf{r}) = -\mathbf{n} \delta_s, \qquad (5) $$

with $\delta_s$ the Dirac delta-function.

For any point inside the volume $V$, bounding surface excluded, Eq. (4) restitutes with $f(\mathbf{r}) = 1$ and $\nabla f(\mathbf{r}) = 0$ the usual equation of motion given in Eq. (1). To see how the $f(\mathbf{r})$ function ensures the stress-free boundary condition on the outer surface, it is necessary to come back to the divergence of stress. Let us evaluate it on an elementary box of volume $V_{\text{box}}$ as shown in Fig. 1.

The top surface of the box $S_t$ belongs to the stress-free surface and the box thickness is infinitesimal. $S_s$ and $S_b$ denote respectively the side and bottom surfaces. Let us evaluate Eq. (4) with integration over the volume of the box

$$ \iiint_{V_{\text{box}}} \{[\nabla f(\mathbf{r})] \cdot \mathbf{T}(\mathbf{r},t) + f(\mathbf{r}) [\nabla \cdot \mathbf{T}(\mathbf{r},t)]\} dV = \iiint_{V_{\text{box}}} \rho(\mathbf{r}) f(\mathbf{r}) \frac{\partial^2 \mathbf{u}(\mathbf{r},t)}{\partial t^2} dV. \quad (6) $$

The first term of the left-hand side, by taking advantage of $\delta_s$, gives

$$ \iiint_{V_{\text{box}}} [\nabla f(\mathbf{r})] \cdot \mathbf{T}(\mathbf{r},t) dV = \iiint_{V_{\text{box}}} -\mathbf{n} \delta_s \cdot \mathbf{T}(\mathbf{r},t) dV = \iint_{S_t} -\mathbf{n} \mathbf{T}(\mathbf{r},t) dS $$

and the second term by applying the divergence theorem

$$ \iiint_{V_{\text{box}}} f(\mathbf{r}) [\nabla \cdot \mathbf{T}(\mathbf{r},t)] dV = \iint_{S_t} \mathbf{n} \mathbf{T}(\mathbf{r},t) dS + \iint_{S_b} -\mathbf{n} \mathbf{T}(\mathbf{r},t) dS + \iint_{S_s} \mathbf{T}(\mathbf{r},t) dS. $$

The contribution from the side of the box is negligible since its height is infinitesimal. The right-hand side of Eq. (6) would be zero since the volume of the box is infinitely small. Thus Eq. (6) leads to

$$ \iint_{S_b} -\mathbf{n} \mathbf{T}(\mathbf{r},t) dS = 0. \qquad (7) $$

Assuming that the cross-sectional area of the box is small enough that the stress is constant across the bottom surface then the normal stress vanishes on and near the bounding surface.

In conclusion, thanks to the function $f(\mathbf{r})$, Eq. (3) incorporates both Newton's second law Eq. (1) and the stress-free boundary condition.

## B. Electrically free boundary condition

For piezoelectric materials, the same device can be applied with the Maxwell equation

$$ \nabla \cdot \mathbf{D} = 0 \qquad (8) $$

to impose, when the electric fields outside the body are assumed to vanish, the electric boundary condition $D_n = 0$ for the normal electric displacement $D_n$ at the electrically free bounding surface.

In theory the electric fields are nonzero outside the body but the error introduced by this assumption is small because most piezoelectric materials have a high dielectric constant.²² For

<figure><img src="image_9.png" /><figcaption>FIG. 1. Elementary box for the evaluation of the divergence of stress.</figcaption></figure>
simple geometries, for instance semi-infinite crystals or plates, it happens that the analytical solving of Laplace's equation is possible and the normal electric displacement in free-space, $D_n$ (free space), can be calculated. It is then easy to avoid the above approximation $D_n = 0$ by using $D'_n = [D_n - D_n(\text{free space})]$ rather than $D_n$ in the Maxwell Eq. (8). This has been largely used with matrix methods³⁶ to impose an exact electrical boundary condition each time the free space electric displacement can be calculated by the Laplace's equation.

In conclusion, here also thanks to the function $f(\mathbf{r})$, Gauss' law Eq. (8) incorporates the electric boundary condition on the bounding surface.

### C. Externally applied normal mechanical force

The case of an externally applied normal mechanical force can be easily taken into account. If we replace the first term on the left-hand side of Eq. (4) $[\nabla f(\mathbf{r})] \cdot T(\mathbf{r}, t)$ by $[\nabla f(\mathbf{r})] \cdot \{T(\mathbf{r}, t) - T_{\text{ext}}(\mathbf{r}, t)\}$ where $T_{\text{ext}}(\mathbf{r}, t)$ denotes the external mechanical boundary excitation, Eq. (6) leads to $T(\mathbf{r}, t) = T_{\text{ext}}(\mathbf{r}, t)$ on and near the bounding surface, which reveals the contribution of the external mechanical source.

### D. Position-dependent physical constants

The function $f(\mathbf{r})$ can be incorporated into the motion equations by means of position-dependent physical constants in the constitutive relations,³⁷ stiffness tensor $c$, piezoelectric tensor $e$, permittivity tensor $\epsilon$, and mass density $\rho$

$$T_i(\mathbf{r}) = c_{ij}^{\text{str}}(\mathbf{r}) S_j(\mathbf{r}) - e_{ij}^{\text{str}}(\mathbf{r}) E_j(\mathbf{r}), \qquad (9a)$$

$$D_i(\mathbf{r}) = \epsilon_{ij}^{\text{str}}(\mathbf{r}) E_j(\mathbf{r}) + e_{ij}^{\text{str}}(\mathbf{r}) S_j(\mathbf{r}), \qquad (9b)$$

with $S_j(\mathbf{r})$ the components of strain tensor and $E_j(\mathbf{r})$ the components of electric field. The superscript $\text{str}$ represents the effective quantity of the structure.³⁸ The contracted Voigt notation is used. Summation over repeated indices is implied throughout this paper. It is assumed that the position-dependent physical constants are nonvanishing only inside the body.

If we consider a uniform material body with a stress-free bounding surface, this gives

$$c_{ij}^{\text{str}}(\mathbf{r}) = c_{ij} \cdot f(\mathbf{r}),$$

$$e_{ij}^{\text{str}}(\mathbf{r}) = e_{ij} \cdot f(\mathbf{r}),$$

$$\epsilon_{ij}^{\text{str}}(\mathbf{r}) = \epsilon_{ij} \cdot f(\mathbf{r}),$$

$$\rho^{\text{str}}(\mathbf{r}) = \rho \cdot f(\mathbf{r}),$$

where $f(\mathbf{r})$ describes the region occupied by the body and $c_{ij}$, $e_{ij}$, $\epsilon_{ij}$ and $\rho$ are the physical constants of the uniform material constituting the body. The constitutive equations are then written as

$$T_i(\mathbf{r}) = [c_{ij}S_j(\mathbf{r}) - e_{ij}E_j(\mathbf{r})] f(\mathbf{r}), \qquad (10a)$$

$$D_i(\mathbf{r}) = [\epsilon_{ij}E_j(\mathbf{r}) + e_{ij}S_j(\mathbf{r})] f(\mathbf{r}). \qquad (10b)$$

For a non uniform body, multilayered or graded for instance, the device of position-dependent physical constants¹⁷,³⁰,³⁸ can also be easily dealt with and the constitutive equations become

$$T_i(\mathbf{r}) = [c_{ij}(\mathbf{r}) S_j(\mathbf{r}) - e_{ij}(\mathbf{r}) E_j(\mathbf{r})] f(\mathbf{r}),$$

$$D_i(\mathbf{r}) = [\epsilon_{ij}(\mathbf{r}) E_j(\mathbf{r}) + e_{ij}(\mathbf{r}) S_j(\mathbf{r})] f(\mathbf{r}),$$

where $c_{ij}(\mathbf{r})$, $e_{ij}(\mathbf{r})$ and $\epsilon_{ij}(\mathbf{r})$ describe the materials and $f(\mathbf{r})$ the region occupied by the structure.

## III. METHOD OF SOLUTION

The mapped orthogonal functions method works whatever the coordinate system adopted, in general chosen to suit to the geometry of the problem under study. So in the constitutive relations
simple geometries, for instance semi-infinite crystals or plates, it happens that the analytical solving of Laplace's equation is possible and the normal electric displacement in free-space, $D_n$ (free space), can be calculated. It is then easy to avoid the above approximation $D_n = 0$ by using $D'_n = [D_n - D_n(\text{free space})]$ rather than $D_n$ in the Maxwell Eq. (8). This has been largely used with matrix methods³⁶ to impose an exact electrical boundary condition each time the free space electric displacement can be calculated by the Laplace's equation.

In conclusion, here also thanks to the function $f(\mathbf{r})$, Gauss' law Eq. (8) incorporates the electric boundary condition on the bounding surface.

### C. Externally applied normal mechanical force

The case of an externally applied normal mechanical force can be easily taken into account. If we replace the first term on the left-hand side of Eq. (4) $[\nabla f(\mathbf{r})] \cdot T(\mathbf{r},t)$ by $[\nabla f(\mathbf{r})] \cdot \{T(\mathbf{r},t) - T_{\text{ext}}(\mathbf{r},t)\}$ where $T_{\text{ext}}(\mathbf{r},t)$ denotes the external mechanical boundary excitation, Eq. (6) leads to $T(\mathbf{r},t) = T_{\text{ext}}(\mathbf{r},t)$ on and near the bounding surface, which reveals the contribution of the external mechanical source.

### D. Position-dependent physical constants

The function $f(\mathbf{r})$ can be incorporated into the motion equations by means of position-dependent physical constants in the constitutive relations,³⁷ stiffness tensor $c$, piezoelectric tensor $e$, permittivity tensor $\epsilon$, and mass density $\rho$

$$T_i(\mathbf{r}) = c_{ij}^{\text{str}}(\mathbf{r}) S_j(\mathbf{r}) - e_{ij}^{\text{str}}(\mathbf{r}) E_j(\mathbf{r}), \qquad (9a)$$

$$D_i(\mathbf{r}) = \epsilon_{ij}^{\text{str}}(\mathbf{r}) E_j(\mathbf{r}) + e_{ij}^{\text{str}}(\mathbf{r}) S_j(\mathbf{r}), \qquad (9b)$$

with $S_j(\mathbf{r})$ the components of strain tensor and $E_j(\mathbf{r})$ the components of electric field. The superscript $\text{str}$ represents the effective quantity of the structure.³⁸ The contracted Voigt notation is used. Summation over repeated indices is implied throughout this paper. It is assumed that the position-dependent physical constants are nonvanishing only inside the body.

If we consider a uniform material body with a stress-free bounding surface, this gives

$$c_{ij}^{\text{str}}(\mathbf{r}) = c_{ij} \cdot f(\mathbf{r}),$$

$$e_{ij}^{\text{str}}(\mathbf{r}) = e_{ij} \cdot f(\mathbf{r}),$$

$$\epsilon_{ij}^{\text{str}}(\mathbf{r}) = \epsilon_{ij} \cdot f(\mathbf{r}),$$

$$\rho^{\text{str}}(\mathbf{r}) = \rho \cdot f(\mathbf{r}),$$

where $f(\mathbf{r})$ describes the region occupied by the body and $c_{ij}$, $e_{ij}$, $\epsilon_{ij}$ and $\rho$ are the physical constants of the uniform material constituting the body. The constitutive equations are then written as

$$T_i(\mathbf{r}) = [c_{ij}S_j(\mathbf{r}) - e_{ij}E_j(\mathbf{r})] f(\mathbf{r}), \qquad (10a)$$

$$D_i(\mathbf{r}) = [\epsilon_{ij}E_j(\mathbf{r}) + e_{ij}S_j(\mathbf{r})] f(\mathbf{r}). \qquad (10b)$$

For a non uniform body, multilayered or graded for instance, the device of position-dependent physical constants¹⁷,³⁰,³⁸ can also be easily dealt with and the constitutive equations become

$$T_i(\mathbf{r}) = [c_{ij}(\mathbf{r}) S_j(\mathbf{r}) - e_{ij}(\mathbf{r}) E_j(\mathbf{r})] f(\mathbf{r}),$$

$$D_i(\mathbf{r}) = [\epsilon_{ij}(\mathbf{r}) E_j(\mathbf{r}) + e_{ij}(\mathbf{r}) S_j(\mathbf{r})] f(\mathbf{r}),$$

where $c_{ij}(\mathbf{r})$, $e_{ij}(\mathbf{r})$ and $\epsilon_{ij}(\mathbf{r})$ describe the materials and $f(\mathbf{r})$ the region occupied by the structure.

## III. METHOD OF SOLUTION

The mapped orthogonal functions method works whatever the coordinate system adopted, in general chosen to suit to the geometry of the problem under study. So in the constitutive relations
Eqs. (9a) and (9b), the strain components must be expressed with respect to the chosen coordinate system.

In rectangular Cartesian coordinates with $\mathbf{r} = (x,y,z)$ and $\mathbf{u} = (u_x, u_y, u_z)$, the strain components are

$$ \begin{pmatrix} S_1 \\ S_2 \\ S_3 \\ S_4 \\ S_5 \\ S_6 \end{pmatrix} = \begin{pmatrix} \frac{\partial u_x}{\partial x} \\ \frac{\partial u_y}{\partial y} \\ \frac{\partial u_z}{\partial z} \\ \frac{\partial u_y}{\partial z} + \frac{\partial u_z}{\partial y} \\ \frac{\partial u_x}{\partial z} + \frac{\partial u_z}{\partial x} \\ \frac{\partial u_x}{\partial y} + \frac{\partial u_y}{\partial x} \end{pmatrix}, \qquad (11a) $$

in cylindrical coordinates with $\mathbf{r} = (r, \theta, z)$ and $\mathbf{u} = (u_r, u_\theta, u_z)$

$$ \begin{pmatrix} S_1 \\ S_2 \\ S_3 \\ S_4 \\ S_5 \\ S_6 \end{pmatrix} = \begin{pmatrix} \frac{\partial u_r}{\partial r} \\ \frac{u_r}{r} + \frac{1}{r} \frac{\partial u_\theta}{\partial \theta} \\ \frac{\partial u_z}{\partial z} \\ \frac{\partial u_\theta}{\partial z} + \frac{1}{r} \frac{\partial u_z}{\partial \theta} \\ \frac{\partial u_r}{\partial z} + \frac{\partial u_z}{\partial r} \\ \frac{1}{r} \frac{\partial u_r}{\partial \theta} + \frac{\partial u_\theta}{\partial r} - \frac{u_\theta}{r} \end{pmatrix} \qquad (11b) $$

and in spherical coordinates with $\mathbf{r} = (r, \theta, \varphi)$ and $\mathbf{u} = (u_r, u_\theta, u_\varphi)$

$$ \begin{pmatrix} S_1 \\ S_2 \\ S_3 \\ S_4 \\ S_5 \\ S_6 \end{pmatrix} = \begin{pmatrix} \frac{\partial u_r}{\partial r} \\ \frac{u_r}{r} + \frac{1}{r} \frac{\partial u_\theta}{\partial \theta} \\ \frac{u_r}{r} + \frac{\cot(\theta)}{r} u_\theta + \frac{1}{r \sin(\theta)} \frac{\partial u_\varphi}{\partial \varphi} \\ \frac{1}{r \sin(\theta)} \frac{\partial u_\theta}{\partial \varphi} + \frac{1}{r} \frac{\partial u_\varphi}{\partial \theta} - \frac{\cot(\theta)}{r} u_\varphi \\ \frac{1}{r \sin(\theta)} \frac{\partial u_r}{\partial \varphi} + \frac{\partial u_\varphi}{\partial r} - \frac{u_\varphi}{r} \\ \frac{1}{r} \frac{\partial u_r}{\partial \theta} + \frac{\partial u_\theta}{\partial r} - \frac{u_\theta}{r} \end{pmatrix}. \qquad (11c) $$

In what follows, in order to keep generality irrespective of the coordinate system, we denote $r_1, r_2$ and $r_3$ the coordinates.

The studies of acoustic wave-based devices making use of position-dependent physical constants have recourse to a mapped orthogonal functions method developed by A.A. Maradudin and his colleagues.³⁹ We will focus on two types of structures:
Eqs. (9a) and (9b), the strain components must be expressed with respect to the chosen coordinate system.

In rectangular Cartesian coordinates with $\mathbf{r} = (x,y,z)$ and $\mathbf{u} = (u_x, u_y, u_z)$, the strain components are

$$ \begin{pmatrix} S_1 \\ S_2 \\ S_3 \\ S_4 \\ S_5 \\ S_6 \end{pmatrix} = \begin{pmatrix} \frac{\partial u_x}{\partial x} \\ \frac{\partial u_y}{\partial y} \\ \frac{\partial u_z}{\partial z} \\ \frac{\partial u_y}{\partial z} + \frac{\partial u_z}{\partial y} \\ \frac{\partial u_x}{\partial z} + \frac{\partial u_z}{\partial x} \\ \frac{\partial u_x}{\partial y} + \frac{\partial u_y}{\partial x} \end{pmatrix}, \qquad (11a) $$

in cylindrical coordinates with $\mathbf{r} = (r, \theta, z)$ and $\mathbf{u} = (u_r, u_\theta, u_z)$

$$ \begin{pmatrix} S_1 \\ S_2 \\ S_3 \\ S_4 \\ S_5 \\ S_6 \end{pmatrix} = \begin{pmatrix} \frac{\partial u_r}{\partial r} \\ \frac{u_r}{r} + \frac{1}{r} \frac{\partial u_\theta}{\partial \theta} \\ \frac{\partial u_z}{\partial z} \\ \frac{\partial u_\theta}{\partial z} + \frac{1}{r} \frac{\partial u_z}{\partial \theta} \\ \frac{\partial u_r}{\partial z} + \frac{\partial u_z}{\partial r} \\ \frac{1}{r} \frac{\partial u_r}{\partial \theta} + \frac{\partial u_\theta}{\partial r} - \frac{u_\theta}{r} \end{pmatrix} \qquad (11b) $$

and in spherical coordinates with $\mathbf{r} = (r, \theta, \varphi)$ and $\mathbf{u} = (u_r, u_\theta, u_\varphi)$

$$ \begin{pmatrix} S_1 \\ S_2 \\ S_3 \\ S_4 \\ S_5 \\ S_6 \end{pmatrix} = \begin{pmatrix} \frac{\partial u_r}{\partial r} \\ \frac{u_r}{r} + \frac{1}{r} \frac{\partial u_\theta}{\partial \theta} \\ \frac{u_r}{r} + \frac{\cot(\theta)}{r} u_\theta + \frac{1}{r \sin(\theta)} \frac{\partial u_\varphi}{\partial \varphi} \\ \frac{1}{r \sin(\theta)} \frac{\partial u_\theta}{\partial \varphi} + \frac{1}{r} \frac{\partial u_\varphi}{\partial \theta} - \frac{\cot(\theta)}{r} u_\varphi \\ \frac{1}{r \sin(\theta)} \frac{\partial u_r}{\partial \varphi} + \frac{\partial u_\varphi}{\partial r} - \frac{u_\varphi}{r} \\ \frac{1}{r} \frac{\partial u_r}{\partial \theta} + \frac{\partial u_\theta}{\partial r} - \frac{u_\theta}{r} \end{pmatrix}. \qquad (11c) $$

In what follows, in order to keep generality irrespective of the coordinate system, we denote $r_1, r_2$ and $r_3$ the coordinates.

The studies of acoustic wave-based devices making use of position-dependent physical constants have recourse to a mapped orthogonal functions method developed by A.A. Maradudin and his colleagues.³⁹ We will focus on two types of structures:
* first type relating to a structure with propagation of the guided wave in a direction given by the wave vector k. This first type will be illustrated on a multilayered plate,

* and a second type relating to structures with stationary wave. This second type will be illustrated on both a 2D rectangular resonator and a 3D cylindrical resonator with transversal isotropy. We adopt the quasistatic approximation

$$E = -\nabla\phi), \quad (12)$$

where $\phi$ denotes the electric potential.

## A. Acoustic wave-based devices with guided wave, $r_1$ propagation

To obtain the dispersion relation, it is necessary to solve both the equation of motion and the Maxwell equation respectively Eq. (1) and Eq. (8) which govern the wave propagation in the device. For the steady-state solution, we assume the mechanical displacement components $u_a$ and the electric potential $\phi$ are of the form

$$u_a(r_1, r_2, r_3, t) = U_a(r_2, r_3) \exp[j(kr_1 - \omega t)], \quad (13a)$$

$$\phi(r_1, r_2, r_3, t) = \Phi(r_2, r_3) \exp[j(kr_1 - \omega t)], \quad (13b)$$

where $U_a$ and $\Phi$ are the magnitudes of the fields and $\omega$ the radian frequency. The subscript $a$ takes on the values 1, 2 and 3. The wave propagates in the $r_1$ direction and the structure is assumed infinite in that direction. For the other two directions $r_2$ and $r_3$ we assume $r_{2\min} \le r_2 \le r_{2\max}$ and $r_{3\min} \le r_3 \le r_{3\max}$ without restrictions on the values of $r_{2\min}$, $r_{2\max}$, $r_{3\min}$ and $r_{3\max}$, some of which can be infinite.

The field quantities are assumed to exist only inside the body so they have to be expanded in a double series of polynomials $P_m(r_2)$ and $Q_n(r_3)$ chosen orthonormal and complete in the respective regions $r_{2\min} \le r_2 \le r_{2\max}$ and $r_{3\min} \le r_3 \le r_{3\max}$

$$U_a(r_2, r_3) = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} p_{mn}^a P_m(r_2) Q_n(r_3), \quad (14a)$$

$$\Phi(r_2, r_3) = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} r_{mn} P_m(r_2) Q_n(r_3), \quad (14b)$$

where $p_{mn}^a$ and $r_{mn}$ are the expansion coefficients. Even if not compulsory and not presented here in the interest of simplification, it is interesting to map the space intervals $[r_{2\min}, r_{2\max}]$ and $[r_{3\min}, r_{3\max}]$ to the orthogonality intervals of the classical orthogonal polynomials, for instance $[-1, +1]$ for Legendre polynomials and $[0, +\infty[$ for Laguerre polynomials. The steps of the solving are as follows:

i) Eqs. (14a) and (14b) are substituted:

a. In the relevant strain-displacement relation Eq. (11a) or (11b) or (11c) depending on the involved coordinate system, this gives the strain components.

b. In Eq. (12) which gives the electric field components.

ii) Strain and electric fields components are substituted in the constitutive relations Eqs. (9a) and (9b). The position-dependent physical constants $c_{ij}(r_2, r_3)$, $e_{ij}(r_2, r_3)$ and $\varepsilon_{ij}(r_2, r_3)$ along with the function $f(r)$ describe both the materials and the geometry of the acoustic wave-based device. For illustration purposes, for a homogeneous medium with a stress-free open-circuit bounding surface occupying the infinite half space given by the equation $r_2 \ge 0$, $f(r) = H(r_2)$ where H denotes the unit step function and for a homogeneous medium under the form of a stress-free open-circuit plate occupying the space $0 \le r_2 \le h$, $f(r) = \pi_{0,h}(r_2)$ where $\pi_{0,h}(r_2)$ denotes the window rectangular function defined as $\pi_{0,h}(r_2) = \begin{cases} 1 & \text{if } 0 \le r_2 \le h \\ 0 & \text{elsewhere} \end{cases}$. In case of a short-circuit bounding surface, there is no constraint on the
* first type relating to a structure with propagation of the guided wave in a direction given by the wave vector k. This first type will be illustrated on a multilayered plate,

* and a second type relating to structures with stationary wave. This second type will be illustrated on both a 2D rectangular resonator and a 3D cylindrical resonator with transversal isotropy. We adopt the quasistatic approximation

$$E = -\nabla\phi), \quad (12)$$

where $\phi$ denotes the electric potential.

## A. Acoustic wave-based devices with guided wave, $r_1$ propagation

To obtain the dispersion relation, it is necessary to solve both the equation of motion and the Maxwell equation respectively Eq. (1) and Eq. (8) which govern the wave propagation in the device. For the steady-state solution, we assume the mechanical displacement components $u_a$ and the electric potential $\phi$ are of the form

$$u_a(r_1, r_2, r_3, t) = U_a(r_2, r_3) \exp[j(kr_1 - \omega t)], \quad (13a)$$

$$\phi(r_1, r_2, r_3, t) = \Phi(r_2, r_3) \exp[j(kr_1 - \omega t)], \quad (13b)$$

where $U_a$ and $\Phi$ are the magnitudes of the fields and $\omega$ the radian frequency. The subscript $a$ takes on the values 1, 2 and 3. The wave propagates in the $r_1$ direction and the structure is assumed infinite in that direction. For the other two directions $r_2$ and $r_3$ we assume $r_{2\min} \le r_2 \le r_{2\max}$ and $r_{3\min} \le r_3 \le r_{3\max}$ without restrictions on the values of $r_{2\min}$, $r_{2\max}$, $r_{3\min}$ and $r_{3\max}$, some of which can be infinite.

The field quantities are assumed to exist only inside the body so they have to be expanded in a double series of polynomials $P_m(r_2)$ and $Q_n(r_3)$ chosen orthonormal and complete in the respective regions $r_{2\min} \le r_2 \le r_{2\max}$ and $r_{3\min} \le r_3 \le r_{3\max}$

$$U_a(r_2, r_3) = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} p_{mn}^a P_m(r_2) Q_n(r_3), \quad (14a)$$

$$\Phi(r_2, r_3) = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} r_{mn} P_m(r_2) Q_n(r_3), \quad (14b)$$

where $p_{mn}^a$ and $r_{mn}$ are the expansion coefficients. Even if not compulsory and not presented here in the interest of simplification, it is interesting to map the space intervals $[r_{2\min}, r_{2\max}]$ and $[r_{3\min}, r_{3\max}]$ to the orthogonality intervals of the classical orthogonal polynomials, for instance $[-1, +1]$ for Legendre polynomials and $[0, +\infty[$ for Laguerre polynomials. The steps of the solving are as follows:

i) Eqs. (14a) and (14b) are substituted:

a. In the relevant strain-displacement relation Eq. (11a) or (11b) or (11c) depending on the involved coordinate system, this gives the strain components.

b. In Eq. (12) which gives the electric field components.

ii) Strain and electric fields components are substituted in the constitutive relations Eqs. (9a) and (9b). The position-dependent physical constants $c_{ij}(r_2, r_3)$, $e_{ij}(r_2, r_3)$ and $\varepsilon_{ij}(r_2, r_3)$ along with the function $f(\mathbf{r})$ describe both the materials and the geometry of the acoustic wave-based device. For illustration purposes, for a homogeneous medium with a stress-free open-circuit bounding surface occupying the infinite half space given by the equation $r_2 \ge 0$, $f(\mathbf{r}) = H(r_2)$ where H denotes the unit step function and for a homogeneous medium under the form of a stress-free open-circuit plate occupying the space $0 \le r_2 \le h$, $f(\mathbf{r}) = \pi_{0,h}(r_2)$ where $\pi_{0,h}(r_2)$ denotes the window rectangular function defined as $\pi_{0,h}(r_2) = \begin{cases} 1 & \text{if } 0 \le r_2 \le h \\ 0 & \text{elsewhere} \end{cases}$. In case of a short-circuit bounding surface, there is no constraint on the
normal electric displacement field, so in Eq. (10b) we must remove $f(\mathbf{r})$ and in Eq. (10a), we must keep $f(\mathbf{r})$ to impose the stress-free boundary condition.

iii) Stress and electric displacement components are substituted in Newton's second law Eq. (1) and in Gauss' law Eq. (8).

iv) We then use the orthonormality of the {$P_m(r_2)$} and {$Q_n(r_3)$}: we multiply Newton's second law and Gauss' law by {$P_s(r_2)$} and {$Q_t(r_3)$} and we integrate over $r_2$ from $r_{2\min}$ to $r_{2\max}$ and $r_3$ from $r_{3\min}$ to $r_{3\max}$. This gives rise to the following equations for the expansion coefficients $p_{mn}^a$ and $r_{mn}$

$$ \mathrm{II}_{mnst}^{ab}(k) p_{mn}^{b} + \mathrm{JJ}_{mnst}^{a}(k) r_{mn} = -\rho \omega^2 \mathrm{MM}_{mnst} p_{mn}^{a}, \quad (15a) $$

$$ \mathrm{GG}_{mnst}^{b}(k) p_{mn}^{b} + \mathrm{KK}_{mnst}(k) r_{mn} = 0, \quad (15b) $$

where the matrices $\mathrm{II}$, $\mathrm{JJ}$, $\mathrm{GG}$ and $\mathrm{KK}$ are functions of the wavenumber $k$. Eq. (15b) gives

$$ r_{mn} = -\mathrm{KK}_{mnst}^{-1}(k) \mathrm{GG}_{mnst}^{b}(k) p_{mn}^{b}. \quad (16a) $$

By substituting $r_{mn}$ into Eq. (15a) we get an eigenvalue matrix for the expansion coefficients $p_{mn}^b$

$$ [\mathrm{II}_{mnst}^{ab}(k) - \mathrm{JJ}_{mnst}^{a}(k) \mathrm{KK}_{mnst}^{-1}(k) \mathrm{GG}_{mnst}^{b}(k)] p_{mn}^{b} = -\rho \omega^2 \mathrm{MM}_{mnst} p_{mn}^{a}. \quad (16b) $$

Eq. (16b) can be solved numerically by standard routines with, as input variable, $k$ real and as output variable $\omega$ real. In practice, the infinite summations in Eqs. (14a) and (14b) are truncated at some numbers $M$ and $N$ which yields $0 \le m, s \le M$ and $0 \le s, t \le N$. The dimensionality of the matrix Eq. (16b) is therefore $\{3 \times (M+1) \times (N+1)\} \times \{3 \times (M+1) \times (N+1)\}$. The solutions are assumed to converge when there is no appreciable change in the eigenvalue $\omega$ as $M$ and $N$ are increased.

Dispersion curves $\omega = f(k)$ are easily obtained. The field distributions are calculated as linear combinations of orthonormal polynomials from Eqs. (14a) and (14b), the coefficients $p_{mn}^a$ and $r_{mn}$ being obtained from the eigenvectors in Eqs. (16a) and (16b).

## B. Acoustic wave-based devices with stationary wave

The mapped orthogonal functions method can be used to model devices with stationary wave such as resonators and transformers. We will illustrate this on a 2D rectangular resonator and a 3D cylindrical one. We assume they both consist of a hexagonal uniform piezoelectric layer belonging to the 6 mm class with its c axis along $r_3$. The layer is assumed sandwiched between two thin ideal metal electrodes without stiffness and mass. The voltage applied to the resonator between the upper and lower metallized bases is

$$ V(t) = V_0 e^{j \omega t}. \quad (17) $$

Two analyses can be carried out for stationary guided acoustic wave-based devices: a free vibration analysis which yields resonant frequencies and mode shapes (modal analysis) and a forced vibration analysis which yields electrical response. Both analyses will be addressed.

### 1. 2D rectangular resonator and boundary conditions

The geometry of the 2D rectangular resonator calls for a Cartesian coordinate system ($r_1, r_2, r_3$). We limit our study to a uniform 2D rectangular resonator and we assume a structure infinite in the $r_2$ direction with no variation in that direction. The resonator is defined by the equations $0 \le r_1 \le L$ and $-\frac{h}{2} \le r_3 \le +\frac{h}{2}$ as shown in Fig. 2.

All the surfaces of the rectangular resonator are stress-free so in Eq. (10a)

$$ f(\mathbf{r}) = \pi_{0,L}(r_1) \times \pi_{-\frac{h}{2},+\frac{h}{2}}(r_3). $$
normal electric displacement field, so in Eq. (10b) we must remove $f(\mathbf{r})$ and in Eq. (10a), we must keep $f(\mathbf{r})$ to impose the stress-free boundary condition.

iii) Stress and electric displacement components are substituted in Newton's second law Eq. (1) and in Gauss' law Eq. (8).

iv) We then use the orthonormality of the {$P_m(r_2)$} and {$Q_n(r_3)$}: we multiply Newton's second law and Gauss' law by {$P_s(r_2)$} and {$Q_t(r_3)$} and we integrate over $r_2$ from $r_{2\min}$ to $r_{2\max}$ and $r_3$ from $r_{3\min}$ to $r_{3\max}$. This gives rise to the following equations for the expansion coefficients $p_{mn}^a$ and $r_{mn}$

$$II_{mnst}^{ab}(k) p_{mn}^{b} + JJ_{mnst}^{a}(k) r_{mn} = -\rho \omega^2 MM_{mnst} p_{mn}^{a}, \quad (15a)$$

$$GG_{mnst}^{b}(k) p_{mn}^{b} + KK_{mnst}(k) r_{mn} = 0, \quad (15b)$$

where the matrices $II$, $JJ$, $GG$ and $KK$ are functions of the wavenumber $k$. Eq. (15b) gives

$$r_{mn} = -KK_{mnst}^{-1}(k) GG_{mnst}^{b}(k) p_{mn}^{b}. \quad (16a)$$

By substituting $r_{mn}$ into Eq. (15a) we get an eigenvalue matrix for the expansion coefficients $p_{mn}^b$

$$[II_{mnst}^{ab}(k) - JJ_{mnst}^{a}(k) KK_{mnst}^{-1}(k) GG_{mnst}^{b}(k)] p_{mn}^{b} = -\rho \omega^2 MM_{mnst} p_{mn}^{a}. \quad (16b)$$

Eq. (16b) can be solved numerically by standard routines with, as input variable, $k$ real and as output variable $\omega$ real. In practice, the infinite summations in Eqs. (14a) and (14b) are truncated at some numbers $M$ and $N$ which yields $0 \le m, s \le M$ and $0 \le s, t \le N$. The dimensionality of the matrix Eq. (16b) is therefore $\{3 \times (M+1) \times (N+1)\} \times \{3 \times (M+1) \times (N+1)\}$. The solutions are assumed to converge when there is no appreciable change in the eigenvalue $\omega$ as $M$ and $N$ are increased.

Dispersion curves $\omega = f(k)$ are easily obtained. The field distributions are calculated as linear combinations of orthonormal polynomials from Eqs. (14a) and (14b), the coefficients $p_{mn}^a$ and $r_{mn}$ being obtained from the eigenvectors in Eqs. (16a) and (16b).

## B. Acoustic wave-based devices with stationary wave

The mapped orthogonal functions method can be used to model devices with stationary wave such as resonators and transformers. We will illustrate this on a 2D rectangular resonator and a 3D cylindrical one. We assume they both consist of a hexagonal uniform piezoelectric layer belonging to the 6 mm class with its c axis along $r_3$. The layer is assumed sandwiched between two thin ideal metal electrodes without stiffness and mass. The voltage applied to the resonator between the upper and lower metallized bases is

$$V(t) = V_0 e^{j \omega t}. \quad (17)$$

Two analyses can be carried out for stationary guided acoustic wave-based devices: a free vibration analysis which yields resonant frequencies and mode shapes (modal analysis) and a forced vibration analysis which yields electrical response. Both analyses will be addressed.

### 1. 2D rectangular resonator and boundary conditions

The geometry of the 2D rectangular resonator calls for a Cartesian coordinate system ($r_1, r_2, r_3$). We limit our study to a uniform 2D rectangular resonator and we assume a structure infinite in the $r_2$ direction with no variation in that direction. The resonator is defined by the equations $0 \le r_1 \le L$ and $-\frac{h}{2} \le r_3 \le +\frac{h}{2}$ as shown in Fig. 2.

All the surfaces of the rectangular resonator are stress-free so in Eq. (10a)

$$f(\mathbf{r}) = \pi_{0,L}(r_1) \times \pi_{-\frac{h}{2},+\frac{h}{2}}(r_3).$$
<figure><img src="image_9.png" /><figcaption>FIG. 2. Configuration of the rectangular resonator.</figcaption></figure>

The normal electric displacement field must vanish only on the lateral surfaces so in Eq. (10b)

$$f(\mathbf{r}) = \pi_{0,L}(r_1).$$

## 2. Cylindrical resonator and boundary conditions

We also model a uniform 3D cylindrical resonator shown in Fig. 3. The geometry of the resonator calls this time for a cylindrical coordinate system $(r_1, r_2, r_3)$. The resonator is defined by the equations $0 \le r_1 \le R$ and $-\frac{h}{2} \le r_3 \le +\frac{h}{2}$.

All the surfaces of the cylindrical resonator are stress-free so in Eq. (10a)

$$f(\mathbf{r}) = \Gamma_R(r_1) \times \pi_{-\frac{h}{2},+\frac{h}{2}}(r_3),$$

with

$$\Gamma_R(r_1) = \begin{cases} 1 & \text{if } 0 \le r_1 \le R \\ 0 & \text{elsewhere} \end{cases}.$$

The normal electric displacement field must vanish only on the lateral surface so in Eq. (10b)

$$f(\mathbf{r}) = \Gamma_R(r_1).$$

<figure><img src="image_10.png" /><figcaption>FIG. 3. Configuration of the cylindrical resonator.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 2. Configuration of the rectangular resonator.</figcaption></figure>

The normal electric displacement field must vanish only on the lateral surfaces so in Eq. (10b)

$$f(\mathbf{r}) = \pi_{0,L}(r_1).$$

## 2. Cylindrical resonator and boundary conditions

We also model a uniform 3D cylindrical resonator shown in Fig. 3. The geometry of the resonator calls this time for a cylindrical coordinate system $(r_1, r_2, r_3)$. The resonator is defined by the equations $0 \le r_1 \le R$ and $-\frac{h}{2} \le r_3 \le +\frac{h}{2}$.

All the surfaces of the cylindrical resonator are stress-free so in Eq. (10a)

$$f(\mathbf{r}) = \Gamma_R(r_1) \times \pi_{-\frac{h}{2},+\frac{h}{2}}(r_3),$$

with

$$\Gamma_R(r_1) = \begin{cases} 1 & \text{if } 0 \le r_1 \le R \\ 0 & \text{elsewhere} \end{cases}.$$

The normal electric displacement field must vanish only on the lateral surface so in Eq. (10b)

$$f(\mathbf{r}) = \Gamma_R(r_1).$$

<figure><img src="image_11.png" /><figcaption>FIG. 3. Configuration of the cylindrical resonator.</figcaption></figure>
### 3. Mechanical displacement and electric potential

For both resonators, the mechanical displacement components and the electric potential are chosen of the form

$$u_a(r_1, r_3, t) = U_a(r_1, r_3) e^{j \omega t}, \qquad (18a)$$

$$\Phi(r_1, r_3, t) = \left[ V_0 \frac{r_3}{h} + \left( r_3^2 - \frac{h^2}{4} \right) \Phi(r_1, r_3) \right] e^{j \omega t}. \qquad (18b)$$

The analytical form of the electric potential is chosen in order to incorporate in the equations of motion the electric source $V(t)$ of Eq. (17) applied between the electrodes of the resonator

$$\Phi(r_1, r_3 = +\frac{h}{2}, t) - \Phi(r_1, r_3 = -\frac{h}{2}, t) = V(t).$$

The field quantities are assumed to exist only inside the resonator and are expanded in a double series of orthonormal polynomials in $r_1$ and $r_3$

$$U_a(r_1, r_3) = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} p_{mn}^a P_m(r_1) Q_n(r_3), \qquad (19a)$$

$$\Phi(r_1, r_3) = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} r_{mn} P_m(r_1) Q_n(r_3), \qquad (19b)$$

with $0 \le r_1 \le L$ for the rectangular resonator, $0 \le r_1 \le R$ for the cylindrical resonator and $-\frac{h}{2} \le r_3 \le +\frac{h}{2}$ for both resonators. The sets of polynomials $\{P_m(r_1)\}$ and $\{Q_n(r_3)\}$ must be chosen complete in their respective intervals of the coordinates $r_1$ and $r_3$. The potential $\Phi(r_1, r_3)$ from Eq. (19b), when incorporated into Eq. (18b) will allow, through the $r_{mn}$'s, to fit the electric potential gradient within the resonator.

### 4. Analytical solutions

The steps of the solving are similar to those described for the acoustic-wave-based devices with guided waves except that there is no wave number and no unknown radian frequency, there is only one system solution. Eqs. (15a) and (15b) are replaced by

$$I_{mnst}^{ab} p_{mn}^{b} + J J_{mnst}^{a} r_{mn} + F_{st}^{a} V_{0} = -\rho \omega^{2} M M_{mnst} p_{mn}^{a}, \qquad (20a)$$

$$GG_{mnst}^{b} p_{mn}^{b} + KK_{mnst} r_{mn} + G_{st} V_{0} = 0. \qquad (20b)$$

Substituting Eq. (20b) into Eq. (20a) yields

$$[I_{mnst}^{ab} - J J_{mnst}^{a} K K_{mnst}^{-1} G G_{mnst}^{b}] p_{mn}^{b} + [F_{st}^{a} - J J_{mnst}^{a} K K_{mnst}^{-1} G_{st}] V_{0} = -\rho \omega^{2} M M_{mnst} p_{mn}^{a}. \qquad (21)$$

According to Ampere's law, the displacement current density through the piezoelectric resonator is

$$J = \frac{\partial D_3}{\partial t} = j \omega D_3$$

and the electric current $I_0$ in the resonator of cross-section $S$, spatially averaged through the thickness $h$ of the resonator, is given by

$$I_0 = -j \omega \frac{1}{h} \int_{-\frac{h}{2}}^{+\frac{h}{2}} \iint_S D_3(r_1, r_3) dS dr_3.$$

The electric current is independent of the $r_3$ coordinate. Averaging over $r_3$ improves the numerical convergence process.40 Using the expression of $D_3$ defined in Eq. (9b), we get

$$I_0 = j \omega C (V_0 + G G_{mn}^{b} p_{mn}^{b}), \qquad (22)$$
### 3. Mechanical displacement and electric potential

For both resonators, the mechanical displacement components and the electric potential are chosen of the form

$$u_a(r_1, r_3, t) = U_a(r_1, r_3) e^{j \omega t}, \qquad (18a)$$

$$\Phi(r_1, r_3, t) = \left[ V_0 \frac{r_3}{h} + \left( r_3^2 - \frac{h^2}{4} \right) \Phi(r_1, r_3) \right] e^{j \omega t}. \qquad (18b)$$

The analytical form of the electric potential is chosen in order to incorporate in the equations of motion the electric source V(t) of Eq. (17) applied between the electrodes of the resonator

$$\Phi(r_1, r_3 = +\frac{h}{2}, t) - \Phi(r_1, r_3 = -\frac{h}{2}, t) = V(t).$$

The field quantities are assumed to exist only inside the resonator and are expanded in a double series of orthonormal polynomials in $r_1$ and $r_3$

$$U_a(r_1, r_3) = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} p_{mn}^a P_m(r_1) Q_n(r_3), \qquad (19a)$$

$$\Phi(r_1, r_3) = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} r_{mn} P_m(r_1) Q_n(r_3), \qquad (19b)$$

with $0 \le r_1 \le L$ for the rectangular resonator, $0 \le r_1 \le R$ for the cylindrical resonator and $-\frac{h}{2} \le r_3 \le +\frac{h}{2}$ for both resonators. The sets of polynomials $\{P_m(r_1)\}$ and $\{Q_n(r_3)\}$ must be chosen complete in their respective intervals of the coordinates $r_1$ and $r_3$. The potential $\Phi(r_1, r_3)$ from Eq. (19b), when incorporated into Eq. (18b) will allow, through the $r_{mn}$'s, to fit the electric potential gradient within the resonator.

### 4. Analytical solutions

The steps of the solving are similar to those described for the acoustic-wave-based devices with guided waves except that there is no wave number and no unknown radian frequency, there is only one system solution. Eqs. (15a) and (15b) are replaced by

$$I_{mnst}^{ab} p_{mn}^{b} + J J_{mnst}^{a} r_{mn} + F_{st}^{a} V_{0} = - \rho \omega^{2} M M_{mnst} p_{mn}^{a}, \qquad (20a)$$

$$GG_{mnst}^{b} p_{mn}^{b} + KK_{mnst} r_{mn} + G_{st} V_{0} = 0. \qquad (20b)$$

Substituting Eq. (20b) into Eq. (20a) yields

$$[I_{mnst}^{ab} - J J_{mnst}^{a} K K_{mnst}^{-1} G G_{mnst}^{b}] p_{mn}^{b} + [F_{st}^{a} - J J_{mnst}^{a} K K_{mnst}^{-1} G_{st}] V_{0} = - \rho \omega^{2} M M_{mnst} p_{mn}^{a}. \qquad (21)$$

According to Ampere's law, the displacement current density through the piezoelectric resonator is

$$J = \frac{\partial D_3}{\partial t} = j \omega D_3$$

and the electric current $I_0$ in the resonator of cross-section S, spatially averaged through the thickness h of the resonator, is given by

$$I_0 = -j \omega \frac{1}{h} \int_{-\frac{h}{2}}^{+\frac{h}{2}} \iint_S D_3(r_1, r_3) dS dr_3.$$

The electric current is independent of the $r_3$ coordinate. Averaging over $r_3$ improves the numerical convergence process.40 Using the expression of $D_3$ defined in Eq. (9b), we get

$$I_0 = j \omega C (V_0 + G G_{mn}^{b} p_{mn}^{b}), \qquad (22)$$
with C the clamped capacitance of the resonator

$$C = \frac{\epsilon_{33}S}{h}.$$

a. Harmonic analysis. The objective is to calculate the electric input admittance of the resonator. Using Eq. (22), we get

$$Y = \frac{I_0}{V_0} = j \omega C \left[ 1 + \frac{1}{V_0} GG_{mn}^b p_{mn}^b \right]. \qquad (23)$$

b. Modal analysis. Our goal here is to calculate the resonance and anti-resonance frequencies and the associated field profiles. The resonance frequencies $\omega_r$ are obtained by vanishing the voltage source, $V_0 = 0$, in Eq. (21)

$$[II_{mnst}^{ab} - JJ_{mnst}^{a} KK_{mnst}^{-1} GG_{mnst}^{b}] p_{mn}^{b} = -\rho \omega_{r}^{2} MM_{mnst} p_{mn}^{a}.$$

The antiresonance frequencies $\omega_a$ are obtained by vanishing the current $I_0 = 0$.

Eq. (22) gives

$$V_0 = -GG_{mn}^b p_{mn}^b. \qquad (24)$$

Substituting Eq. (24) into Eq. (21) yields

$$[II_{mnst}^{ab} - JJ_{mnst}^{a} KK_{mnst}^{-1} GG_{mnst}^{b} - F_{st}^{a} GG_{mn}^{b} + JJ_{mnst}^{a} KK_{mnst}^{-1} G_{st} GG_{mn}^{b}] p_{mn}^{b} = -\rho \omega_{a}^{2} MM_{mnst} p_{mn}^{a}.$$

Here again, in practice the infinite summations in Eqs. (19a) and (19b) are truncated with $0 \le m, s \le M$ and $0 \le s, t \le N$. In both cases, the eigenvectors $p_{mn}^b$ allow the field profiles associated to the resonance and antiresonance frequencies to be calculated.

## IV. ILLUSTRATIONS

The illustrations do not aim at validating the mapped orthogonal functions method as this has already been done through several literature papers for both acoustic wave-based devices with guided waves$^{34,38,39,41}$ and with stationary waves.$^{30,31,42}$ The objective is rather to illustrate in a same paper the potentialities of the method for various geometries and various regimes: first a multilayered plate where the method gives, through a single calculation, all types of free ultrasonic waves, Lamb-like and SH waves, second a 2D rectangular resonator which highlights the capability of the mapped orthogonal functions method to give the forced response which includes both the thickness and lateral parasitic resonances, the latter so important in a design process, and third a 3D cylindrical resonator.

The physical parameters for the materials used for the illustrations are given in Table I.

The mapped orthogonal functions method does not suffer from the large $fd$ problem.

### A. Multilayered plate$^{43}$

The multilayered plate consists of a topmost layer of AlAs followed by three GaAs/AlAs bilayers as shown in Fig. 4.

The effective stiffness tensors are expressed as

$$c_{ij}^{\text{str}}(\mathbf{r}) = \left\{ c_{ij}^{\text{AlAs}} \pi_{0,h_I}(r_3) + \sum_{n=1}^{3} \left[ c_{ij}^{\text{GaAs}} \pi_{h_{2n-1}, h_{2n}}(r_3) + c_{ij}^{\text{AlAs}} \pi_{h_{2n}, h_{2n+1}}(r_3) \right] \right\} \pi_{0,h_7}(r_3).$$

The effective piezoelectric, permittivity tensors and mass density can be expressed similarly. All layers are of equal thickness. The outmost surfaces are open-circuited surfaces. Figs. 5 and 6 present the dispersion curves, phase velocity as a function of the frequency-thickness product $f h$ where $f$ denotes the frequency and $h = h_7$ the plate thickness for respectively the Lamb-like and SH plate modes in the seven-layered plate. The waves are propagating in the [110] direction and the plates surfaces are perpendicular to the Z crystallographic axis [001].
with C the clamped capacitance of the resonator

$$C = \frac{\epsilon_{33}S}{h}.$$

a. Harmonic analysis. The objective is to calculate the electric input admittance of the resonator. Using Eq. (22), we get

$$Y = \frac{I_0}{V_0} = j \omega C \left[ 1 + \frac{1}{V_0} GG_{mn}^b p_{mn}^b \right]. \qquad (23)$$

b. Modal analysis. Our goal here is to calculate the resonance and anti-resonance frequencies and the associated field profiles. The resonance frequencies $\omega_r$ are obtained by vanishing the voltage source, $V_0 = 0$, in Eq. (21)

$$[II_{mnst}^{ab} - JJ_{mnst}^{a} KK_{mnst}^{-1} GG_{mnst}^{b}] p_{mn}^{b} = -\rho \omega_{r}^{2} MM_{mnst} p_{mn}^{a}.$$

The antiresonance frequencies $\omega_a$ are obtained by vanishing the current $I_0 = 0$.

Eq. (22) gives

$$V_0 = -GG_{mn}^b p_{mn}^b. \qquad (24)$$

Substituting Eq. (24) into Eq. (21) yields

$$[II_{mnst}^{ab} - JJ_{mnst}^{a} KK_{mnst}^{-1} GG_{mnst}^{b} - F_{st}^{a} GG_{mn}^{b} + JJ_{mnst}^{a} KK_{mnst}^{-1} G_{st} GG_{mn}^{b}] p_{mn}^{b} = -\rho \omega_{a}^{2} MM_{mnst} p_{mn}^{a}.$$

Here again, in practice the infinite summations in Eqs. (19a) and (19b) are truncated with $0 \le m, s \le M$ and $0 \le s, t \le N$. In both cases, the eigenvectors $p_{mn}^b$ allow the field profiles associated to the resonance and antiresonance frequencies to be calculated.

## IV. ILLUSTRATIONS

The illustrations do not aim at validating the mapped orthogonal functions method as this has already been done through several literature papers for both acoustic wave-based devices with guided waves$^{34,38,39,41}$ and with stationary waves.$^{30,31,42}$ The objective is rather to illustrate in a same paper the potentialities of the method for various geometries and various regimes: first a multilayered plate where the method gives, through a single calculation, all types of free ultrasonic waves, Lamb-like and SH waves, second a 2D rectangular resonator which highlights the capability of the mapped orthogonal functions method to give the forced response which includes both the thickness and lateral parasitic resonances, the latter so important in a design process, and third a 3D cylindrical resonator.

The physical parameters for the materials used for the illustrations are given in Table I.

The mapped orthogonal functions method does not suffer from the large $fd$ problem.

### A. Multilayered plate$^{43}$

The multilayered plate consists of a topmost layer of AlAs followed by three GaAs/AlAs bilayers as shown in Fig. 4.

The effective stiffness tensors are expressed as

$$c_{ij}^{\text{str}}(\mathbf{r}) = \left\{ c_{ij}^{\text{AlAs}} \pi_{0,h_1}(r_3) + \sum_{n=1}^{3} \left[ c_{ij}^{\text{GaAs}} \pi_{h_{2n-1}, h_{2n}(r_3)} + c_{ij}^{\text{AlAs}} \pi_{h_{2n}, h_{2n+1}(r_3)} \right] \right\} \pi_{0,h_7}(r_3).$$

The effective piezoelectric, permittivity tensors and mass density can be expressed similarly. All layers are of equal thickness. The outmost surfaces are open-circuited surfaces. Figs. 5 and 6 present the dispersion curves, phase velocity as a function of the frequency-thickness product $f h$ where $f$ denotes the frequency and $h = h_7$ the plate thickness for respectively the Lamb-like and SH plate modes in the seven-layered plate. The waves are propagating in the [110] direction and the plates surfaces are perpendicular to the Z crystallographic axis [001].
<table><caption>TABLE I. Materials properties.</caption><thead><tr><th>Parameter</th><th>GaAs</th><th>AlAs</th><th>ZnO</th><th>PZT5H</th></tr></thead><tbody><tr><td>ρ (kg m<sup>-3</sup>)</td><td>5307</td><td>3760</td><td>5676</td><td>7500</td></tr><tr><td colspan="5">Stiffness Constants C<sub>ij</sub> (x10<sup>+11</sup> N/m<sup>2</sup>)</td></tr><tr><td>c<sup>E</sup><sub>11</sub></td><td>1.188</td><td>1.202</td><td>2.097</td><td>1.26</td></tr><tr><td>c<sup>E</sup><sub>12</sub></td><td>0.538</td><td>0.570</td><td>1.211</td><td>0.795</td></tr><tr><td>c<sup>E</sup><sub>13</sub></td><td></td><td></td><td>1.051</td><td>0.841</td></tr><tr><td>c<sup>E</sup><sub>33</sub></td><td></td><td></td><td>2.109</td><td>1.17</td></tr><tr><td>c<sup>E</sup><sub>44</sub></td><td>0.594</td><td>0.589</td><td>0.425</td><td>0.23</td></tr><tr><td>c<sup>E</sup><sub>66</sub></td><td></td><td></td><td></td><td>0.2325</td></tr><tr><td colspan="5">Piezoelectric Constants ε<sub>ij</sub> (C.m<sup>-2</sup>)</td></tr><tr><td>e<sub>14</sub></td><td>0.154</td><td>0.225</td><td></td><td></td></tr><tr><td>e<sub>15</sub></td><td></td><td></td><td>-0.59</td><td>17</td></tr><tr><td>e<sub>31</sub></td><td></td><td></td><td>-0.61</td><td>-6.5</td></tr><tr><td>e<sub>33</sub></td><td></td><td></td><td>1.14</td><td>23.3</td></tr><tr><td colspan="5">Relative Permittivity ε<sub>ij</sub>/ε<sub>0</sub></td></tr><tr><td>ε<sup>S</sup><sub>11</sub></td><td>12.5</td><td>10.06</td><td>8.335</td><td>1700</td></tr><tr><td>ε<sup>S</sup><sub>33</sub></td><td></td><td></td><td>8.843</td><td>1470</td></tr></tbody></table>

The structure is symmetric, accordingly the wave motion occurs with either symmetric or antisymmetric field distributions. For the Lamb-like modes, we find two fundamental modes $A_0$ and $S_0$ without cut-off frequency. All higher-order modes end in a cut-off frequency and exhibit a complex behavior where symmetric and antisymmetric modes cross one or several times. For the SH modes, we find one fundamental mode without cut-off frequency while all higher-order modes end in a cut-off frequency.

## B. 2D rectangular resonator

Figure 7 presents as a function of the normalized radian frequency the calculated normalized admittances $Y_{1D}$ and $Y_{2D}$ for two 2µm thick lossy ZnO bulk acoustic wave resonators respectively laterally unbounded and bounded.⁴⁴ The radian frequency is normalized by the fundamental thickness antiresonance radian frequency and the admittance by $jC\omega$ where C denotes the clamped capacitance and $\omega$ the radian frequency. The width of the laterally bounded resonator is 200µm. The losses are incorporated through a complex stiffness component [$c_{33} + j\omega\omega_{33}$] with a viscosity

<figure><img src="image_11.png" /><figcaption>FIG. 4. Description of the multilayered plate.</figcaption></figure>
<table><caption>TABLE I. Materials properties.</caption><thead><tr><th>Parameter</th><th>GaAs</th><th>AlAs</th><th>ZnO</th><th>PZT5H</th></tr></thead><tbody><tr><td>ρ (kg m<sup>-3</sup>)</td><td>5307</td><td>3760</td><td>5676</td><td>7500</td></tr><tr><td colspan="5">Stiffness Constants C<sub>ij</sub> (x10<sup>+11</sup> N/m<sup>2</sup>)</td></tr><tr><td>c<sup>E</sup><sub>11</sub></td><td>1.188</td><td>1.202</td><td>2.097</td><td>1.26</td></tr><tr><td>c<sup>E</sup><sub>12</sub></td><td>0.538</td><td>0.570</td><td>1.211</td><td>0.795</td></tr><tr><td>c<sup>E</sup><sub>13</sub></td><td></td><td></td><td>1.051</td><td>0.841</td></tr><tr><td>c<sup>E</sup><sub>33</sub></td><td></td><td></td><td>2.109</td><td>1.17</td></tr><tr><td>c<sup>E</sup><sub>44</sub></td><td>0.594</td><td>0.589</td><td>0.425</td><td>0.23</td></tr><tr><td>c<sup>E</sup><sub>66</sub></td><td></td><td></td><td></td><td>0.2325</td></tr><tr><td colspan="5">Piezoelectric Constants ε<sub>ij</sub> (C.m<sup>-2</sup>)</td></tr><tr><td>ε<sub>14</sub></td><td>0.154</td><td>0.225</td><td></td><td></td></tr><tr><td>ε<sub>15</sub></td><td></td><td></td><td>-0.59</td><td>17</td></tr><tr><td>ε<sub>31</sub></td><td></td><td></td><td>-0.61</td><td>-6.5</td></tr><tr><td>ε<sub>33</sub></td><td></td><td></td><td>1.14</td><td>23.3</td></tr><tr><td colspan="5">Relative Permittivity ε<sub>ij</sub>/ε<sub>0</sub></td></tr><tr><td>ε<sup>S</sup><sub>11</sub></td><td>12.5</td><td>10.06</td><td>8.335</td><td>1700</td></tr><tr><td>ε<sup>S</sup><sub>33</sub></td><td></td><td></td><td>8.843</td><td>1470</td></tr></tbody></table>

The structure is symmetric, accordingly the wave motion occurs with either symmetric or antisymmetric field distributions. For the Lamb-like modes, we find two fundamental modes $A_0$ and $S_0$ without cut-off frequency. All higher-order modes end in a cut-off frequency and exhibit a complex behavior where symmetric and antisymmetric modes cross one or several times. For the SH modes, we find one fundamental mode without cut-off frequency while all higher-order modes end in a cut-off frequency.

## B. 2D rectangular resonator

Figure 7 presents as a function of the normalized radian frequency the calculated normalized admittances $Y_{1D}$ and $Y_{2D}$ for two 2µm thick lossy ZnO bulk acoustic wave resonators respectively laterally unbounded and bounded.⁴⁴ The radian frequency is normalized by the fundamental thickness antiresonance radian frequency and the admittance by $jC\omega$ where C denotes the clamped capacitance and $\omega$ the radian frequency. The width of the laterally bounded resonator is 200µm. The losses are incorporated through a complex stiffness component [$c_{33} + j\omega\omega_{33}$] with a viscosity

<figure><img src="image_12.png" /><figcaption>FIG. 4. Description of the multilayered plate.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 5. Dispersion curves for Lamb-like plate modes in the seven-layered plate as a function of the frequency-thickness product.</figcaption></figure>

tensor component $\eta_{33} = 11.96 \times 10^{-3}$Pa.s. The mapped orthogonal functions method restitutes for the laterally bounded resonator both the thickness and parasitic lateral resonances.

The mapped orthogonal functions method is also able to retrieve the modal frequencies. Table II gives the first and third resonance and antiresonance frequencies for a 0.61 mm wide and 1.53 mm thick PZT5H resonator calculated by the mapped orthogonal functions method and compared with the experimental and theoretical ones taken from literature.⁴⁵ The agreement is quite good.

<figure><img src="image_13.png" /><figcaption>FIG. 6. Dispersion curves for SH plate modes in the seven-layered plate as a function of the frequency-thickness product.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 5. Dispersion curves for Lamb-like plate modes in the seven-layered plate as a function of the frequency-thickness product.</figcaption></figure>

tensor component $\eta_{33} = 11.96 \times 10^{-3}$Pa.s. The mapped orthogonal functions method restitutes for the laterally bounded resonator both the thickness and parasitic lateral resonances.

The mapped orthogonal functions method is also able to retrieve the modal frequencies. Table II gives the first and third resonance and antiresonance frequencies for a 0.61 mm wide and 1.53 mm thick PZT5H resonator calculated by the mapped orthogonal functions method and compared with the experimental and theoretical ones taken from literature.⁴⁵ The agreement is quite good.

<figure><img src="image_14.png" /><figcaption>FIG. 6. Dispersion curves for SH plate modes in the seven-layered plate as a function of the frequency-thickness product.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIG. 7. Normalized admittance of the rectangular resonator as a function of the normalized radian frequency (M=N=9).</figcaption></figure>

### C. Cylindrical resonator<sup>42</sup>

The normalized electrical input impedance is presented in Fig. 8 for thickness extensional modes and in Fig. 9 for radial modes as a function of the normalized frequency $\Omega$ for a ZnO disc with $D/H = 300$ where $D$ denotes the diameter and $H$ the thickness of the disc. The frequency is normalized by the fundamental thickness resonance frequency. The analytical results are obtained from a one-dimensional structural theory<sup>46</sup> made possible by the transversal isotropy and the high value $D/H = 300$.

The mapped orthogonal functions method is able to retrieve in a single calculation both the thickness and radial modes. They are represented in two separate figures for a good readability of the radial modes.

## V. EXTENSIONS AND POSSIBLE FURTHER DEVELOPMENTS

For the sake of clarity, the method was presented for acoustic wave-based devices which allow to use, for expressing the field quantities, a single set of polynomials over the entire structure. For heterogeneous bodies with regions exhibiting a high contrast, a single set of polynomials gives rise to a slow convergence owing to the fact that the mechanical displacement and the stress are infinitely

<table><caption>TABLE II. Modal frequencies for a PZT5H resonator.</caption><thead><tr><th>Results</th><th>Modes</th><th>Resonance</th><th>Antiresonance</th></tr></thead><tbody><tr><td rowspan="2">Experimental<sup>45</sup></td><td>1</td><td>0.91 MHz</td><td>1.22 MHz</td></tr><tr><td>3</td><td>2.45 MHz</td><td>2.51 MHz</td></tr><tr><td rowspan="2">Theoretical<sup>45</sup></td><td>1</td><td>0.99 MHz</td><td>1.29 MHz</td></tr><tr><td>3</td><td>2.47 MHz</td><td>2.49 MHz</td></tr><tr><td rowspan="2">Mapped orthogonal functions method</td><td>1</td><td>0.94 MHz</td><td>1.26 MHz</td></tr><tr><td>3</td><td>2.48 MHz</td><td>2.52 MHz</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>FIG. 7. Normalized admittance of the rectangular resonator as a function of the normalized radian frequency (M=N=9).</figcaption></figure>

### C. Cylindrical resonator<sup>42</sup>

The normalized electrical input impedance is presented in Fig. 8 for thickness extensional modes and in Fig. 9 for radial modes as a function of the normalized frequency $\Omega$ for a ZnO disc with $D/H = 300$ where $D$ denotes the diameter and $H$ the thickness of the disc. The frequency is normalized by the fundamental thickness resonance frequency. The analytical results are obtained from a one-dimensional structural theory<sup>46</sup> made possible by the transversal isotropy and the high value $D/H = 300$.

The mapped orthogonal functions method is able to retrieve in a single calculation both the thickness and radial modes. They are represented in two separate figures for a good readability of the radial modes.

## V. EXTENSIONS AND POSSIBLE FURTHER DEVELOPMENTS

For the sake of clarity, the method was presented for acoustic wave-based devices which allow to use, for expressing the field quantities, a single set of polynomials over the entire structure. For heterogeneous bodies with regions exhibiting a high contrast, a single set of polynomials gives rise to a slow convergence owing to the fact that the mechanical displacement and the stress are infinitely

<table><caption>TABLE II. Modal frequencies for a PZT5H resonator.</caption><thead><tr><th>Results</th><th>Modes</th><th>Resonance</th><th>Antiresonance</th></tr></thead><tbody><tr><td rowspan="2">Experimental<sup>45</sup></td><td>1</td><td>0.91 MHz</td><td>1.22 MHz</td></tr><tr><td>3</td><td>2.45 MHz</td><td>2.51 MHz</td></tr><tr><td rowspan="2">Theoretical<sup>45</sup></td><td>1</td><td>0.99 MHz</td><td>1.29 MHz</td></tr><tr><td>3</td><td>2.47 MHz</td><td>2.49 MHz</td></tr><tr><td rowspan="2">Mapped orthogonal functions method</td><td>1</td><td>0.94 MHz</td><td>1.26 MHz</td></tr><tr><td>3</td><td>2.48 MHz</td><td>2.52 MHz</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>FIG. 8. Normalized electrical input impedance for thickness extensional modes for a ZnO disc ($D/H = 300, M = N = 10$).</figcaption></figure>

differentiable everywhere in space but at interfaces. This concern has already been addressed previ-
ously.⁹,⁴¹ The high contrast can originate from two contiguous distinct materials with very dissimilar
physical properties or from two regions of the same material with either distinct crystallographic
orientations or one region metallized and the other one non metallized. Adjacent electroded and un-
electroded regions may behave as very dissimilar materials. In fact this issue can be easily overcome
at the cost of a bigger dimensionality of the system matrix. The solution for a correct modelling of
discontinuities, as illustrated in previous papers,⁹,⁴⁷,⁴⁸ consists in using a distinct set of polynomials

<figure><img src="image_16.png" /><figcaption>FIG. 9. Normalized electrical input impedance for radial modes for a ZnO disc ($D/H = 300, M = N = 12$).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIG. 8. Normalized electrical input impedance for thickness extensional modes for a ZnO disc ($D/H = 300, M = N = 10$).</figcaption></figure>

differentiable everywhere in space but at interfaces. This concern has already been addressed previ-
ously.⁹,⁴¹ The high contrast can originate from two contiguous distinct materials with very dissimilar
physical properties or from two regions of the same material with either distinct crystallographic
orientations or one region metallized and the other one non metallized. Adjacent electroded and un-
electroded regions may behave as very dissimilar materials. In fact this issue can be easily overcome
at the cost of a bigger dimensionality of the system matrix. The solution for a correct modelling of
discontinuities, as illustrated in previous papers,⁹,⁴⁷,⁴⁸ consists in using a distinct set of polynomials

<figure><img src="image_17.png" /><figcaption>FIG. 9. Normalized electrical input impedance for radial modes for a ZnO disc ($D/H = 300, M = N = 12$).</figcaption></figure>
for each region of the body with dissimilar materials or for each region with and without electrode
rather than a unique set of polynomials for the entire body. Apart from the study of convergence as a
function of the number of polynomials kept in the series expansions Eqs. (14a) and (14b) for devices
with guided waves and Eqs. (19a) and (19b) for devices with stationary waves, there is no obvious
criterion to decide if a single set of polynomials is appropriate to model a heterogeneous structure.
Indeed this key point would deserve a specific study but the problem is quite difficult. Each material
is characterized by a great number of physical constants, until 45 for a triclinic piezoelectric material
without taking into account other parameters such as mass density, thickness and boundary conditions.
Piezoelectric-piezomagnetic composite structures involve much more physical constants. Based on
the authors' experience, it is important to underline that in a great majority of heterogeneous piezo-
electric structures, a single set of polynomials works quite well except if exact fields discontinuities
at interfaces are needed. We faced difficulties in two specific cases with no convergence: diamond
with silicon or ZnO, diamond always presents a very high contrast with any usual material used in
acoustic wave-based devices, and angle-ply laminated composite plates. Also in partially electroded
resonators, to obtain the parasitic modes, it is recommended to use a distinct set of polynomials for
each region, the electroded one and the non electroded one. In the case of use of more than three
distinct sets of polynomials for a heterogeneous structure, the preparatory mathematical manipulation
to implement the method in a computer code becomes cumbersome and it is recommended to have
recourse to a computer algebra system to facilitate the work. The presented tool is well suited for
modeling film bulk acoustic resonators and much less suitable for solidly mounted resonators which
include a greater number of layers and particularly a high contrast multilayer for the Bragg mirror.

The method was also presented for structures with regular geometries. In fact, nothing prevents
dealing with topographic devices with more complicated shapes. For the resonators, we could for
instance consider curved lateral surfaces defined by an equation of the form r₁ = f (r₃) for some
well-behaved function f (r₃).<sup>17</sup> The main additional task would be the computation of the Dirac and
Heavyside functions with an argument taking into account the relation r₁ = f (r₃).

For cylindrical structures, a further development interesting to develop is for the cases either of
a cylindrical anisotropy more complicated than transverse isotropy or of a rectangular anisotropy:
the mapped orthogonal functions method, with its position-dependent physical constants, allows
to deal simultaneously with the geometry, a complicated anisotropy and eventually a function-
ally graded material. C. Baron claims the same advantage for anisotropy more complicated than
transverse isotropy<sup>13</sup> with the Peano expansion of the matricant based on Stroh's formalism.

Another further development is the calculation of non propagating waves given by complex
solutions. The method such as presented is not able to obtain easily complex solutions and need to
be extended to calculate non propagative waves. The work under progress keeps the advantage of
an eigenvalue system and will be published soon in a subsequent paper. It simply consists in trans-
forming the quadratic eigenvalue problem Eq. (16b) into a linear one via a change of variable<sup>9-11</sup>
with ω real as input variable and k complex as output variable. This results in an increase in the
dimensionality of the matrices by a factor of four for the same number of terms in the expansions
Eqs. (14a) and (14b). Working with a linear eigenvalue system rather than having recourse to a
tedious characteristic equation is crucial when searching for complex solutions.

Moreover, unlike FEM, the mapped orthogonal functions method has not recourse to a mesh
of the structure. So it can be an alternative to FEM for structures size of which calls for too impor-
tant memory resources. For instance, Laguerre polynomials allow to model semi-infinite acoustic
wave-based devices.

A big challenge which remains to face with the mapped orthogonal functions method is to
deal with acoustic wave-based devices including a mechanical or electric source which gives rise to
propagative waves rather than stationary waves that is to say the modelling of acoustic wave-based
radiant devices.

VI. CONCLUSIONS

The mapped orthogonal functions method has been outlined for analyzing various types of
acoustic wave-based devices, some with guided waves and others with stationary waves. The device
for each region of the body with dissimilar materials or for each region with and without electrode
rather than a unique set of polynomials for the entire body. Apart from the study of convergence as a
function of the number of polynomials kept in the series expansions Eqs. (14a) and (14b) for devices
with guided waves and Eqs. (19a) and (19b) for devices with stationary waves, there is no obvious
criterion to decide if a single set of polynomials is appropriate to model a heterogeneous structure.
Indeed this key point would deserve a specific study but the problem is quite difficult. Each material
is characterized by a great number of physical constants, until 45 for a triclinic piezoelectric material
without taking into account other parameters such as mass density, thickness and boundary conditions.
Piezoelectric-piezomagnetic composite structures involve much more physical constants. Based on
the authors' experience, it is important to underline that in a great majority of heterogeneous piezo-
electric structures, a single set of polynomials works quite well except if exact fields discontinuities
at interfaces are needed. We faced difficulties in two specific cases with no convergence: diamond
with silicon or ZnO, diamond always presents a very high contrast with any usual material used in
acoustic wave-based devices, and angle-ply laminated composite plates. Also in partially electroded
resonators, to obtain the parasitic modes, it is recommended to use a distinct set of polynomials for
each region, the electroded one and the non electroded one. In the case of use of more than three
distinct sets of polynomials for a heterogeneous structure, the preparatory mathematical manipulation
to implement the method in a computer code becomes cumbersome and it is recommended to have
recourse to a computer algebra system to facilitate the work. The presented tool is well suited for
modeling film bulk acoustic resonators and much less suitable for solidly mounted resonators which
include a greater number of layers and particularly a high contrast multilayer for the Bragg mirror.

The method was also presented for structures with regular geometries. In fact, nothing prevents
dealing with topographic devices with more complicated shapes. For the resonators, we could for
instance consider curved lateral surfaces defined by an equation of the form r₁ = f (r₃) for some
well-behaved function f (r₃).<sup>17</sup> The main additional task would be the computation of the Dirac and
Heavyside functions with an argument taking into account the relation r₁ = f (r₃).

For cylindrical structures, a further development interesting to develop is for the cases either of
a cylindrical anisotropy more complicated than transverse isotropy or of a rectangular anisotropy:
the mapped orthogonal functions method, with its position-dependent physical constants, allows
to deal simultaneously with the geometry, a complicated anisotropy and eventually a function-
ally graded material. C. Baron claims the same advantage for anisotropy more complicated than
transverse isotropy<sup>13</sup> with the Peano expansion of the matricant based on Stroh's formalism.

Another further development is the calculation of non propagating waves given by complex
solutions. The method such as presented is not able to obtain easily complex solutions and need to
be extended to calculate non propagative waves. The work under progress keeps the advantage of
an eigenvalue system and will be published soon in a subsequent paper. It simply consists in trans-
forming the quadratic eigenvalue problem Eq. (16b) into a linear one via a change of variable<sup>9-11</sup>
with ω real as input variable and k complex as output variable. This results in an increase in the
dimensionality of the matrices by a factor of four for the same number of terms in the expansions
Eqs. (14a) and (14b). Working with a linear eigenvalue system rather than having recourse to a
tedious characteristic equation is crucial when searching for complex solutions.

Moreover, unlike FEM, the mapped orthogonal functions method has not recourse to a mesh
of the structure. So it can be an alternative to FEM for structures size of which calls for too impor-
tant memory resources. For instance, Laguerre polynomials allow to model semi-infinite acoustic
wave-based devices.

A big challenge which remains to face with the mapped orthogonal functions method is to
deal with acoustic wave-based devices including a mechanical or electric source which gives rise to
propagative waves rather than stationary waves that is to say the modelling of acoustic wave-based
radiant devices.

VI. CONCLUSIONS

The mapped orthogonal functions method has been outlined for analyzing various types of
acoustic wave-based devices, some with guided waves and others with stationary waves. The device
to automatically incorporate the mechanical and/or electrical boundary conditions into the equations of motion has been given and explained. Shown are the results for several applications for propagative and stationary waves. Extensions and possible further developments have been addressed.

Position-dependent physical constants are an interesting device which allows to easily incorporate boundary conditions and mechanical and electrical sources into field equations. In addition, they allow to describe exactly heterogeneous and/or inhomogeneous structures.

It is worth mentioning that, for guided waves in a layered plate supported or not by a substrate, the polynomial approach in Ref. 11 and the mapped orthogonal functions method present many similarities. Both use a polynomial expansion, result through an orthogonalization process in a linear eigenvalue problem, can directly determine complex wave numbers for a given frequency and can deal with complicated acoustic wave-based devices, for instance piezoelectric-piezomagnetic composite devices.<fnref n="11" />,⁴⁹ The difference resides in the treatment of the continuity and boundary conditions. For the mapped orthogonal functions method, the continuity of the generalized displacement vector (mechanical displacement components, electric and magnetic potentials) is ensured by a wise choice of the associated polynomial expansion⁴⁷ and the continuity and boundary conditions of the generalized stress vector (normal stress components, normal electric displacement and magnetic induction) are ensured thanks to the δ-functions given by the position-dependent physical constants. All these continuity and boundary conditions are incorporated into equations of motion and the orthogonalization process straightforwardly results in a system with as many equations as unknowns. For (ω real, k complex), the dimensionality of the associated matrix eigenvalue equation for a piezoelectric structure with N<sub>c</sub> homogeneous layers is [2 × N<sub>c</sub> × {3 × (M + 1)}] × [2 × N<sub>c</sub> × {3 × (M + 1)}] where (M + 1) denotes the number of polynomials kept in the expansions of the four components of the generalized displacement vector. For the polynomial approach in Ref. 11 a first part of the system of equations which govern the wave motion is given by the equations imposing the continuity and boundary conditions without orthogonalization process and the second part, through an orthogonalization process, by the equations governing the generalized displacement vector. For (ω real, k complex), the dimensionality of the associated matrix eigenvalue equation for the same piezoelectric structure with N<sub>c</sub> homogeneous layers and the same polynomial expansions is [2 × 4 × N<sub>c</sub> × (M + 1)] × [2 × 4 × N<sub>c</sub> × (M + 1)]. For the eigenvalue problem, the dimensionality is slightly at the advantage of the mapped orthogonal functions method by a factor of nine-sixteenth. The advantage is even more important for piezoelectric-piezomagnetic structures, factor of nine-twenty fifth always at the advantage of the mapped orthogonal functions method. This is counterbalanced, for the approach in Ref. 11, by a lighter preparatory mathematical work when working with more than three distinct sets of polynomials. Probably the polynomial approach in Ref. 11 could be easily extended (i) to deal with through-thickness FGM structures and (ii) to give access, for structures with stationary waves, to modal and forced responses of acoustic wave-based multilayered devices (homogeneous layers supported or not by a substrate).

<footnote n="1">¹ W.T. Thomson, *J. Appl. Phys.* **21**, 89 (1950).</footnote>

<footnote n="2">² L Knopoff, Bull. Seismol. Soc. Am. **54**, 431 (1964).</footnote>

<footnote n="3">³ M.J.S. Lowe, *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **42**(4), 525 (1995).</footnote>

<footnote n="4">⁴ M. Castaings and B. Hosten, in *Ultrasonics International Conference* (Butterworth-Heinemann, Ltd, London, UK, 1993), pp. 431–434.</footnote>

<footnote n="5">⁵ B. Hosten and M. Castaings, *J. Acoust. Soc. Am.* **94**, 1488 (1993).</footnote>

<footnote n="6">⁶ L. Wang and S.I. Rokhlin, *Ultrasonics* **39**, 413 (2001).</footnote>

<footnote n="7">⁷ B. Collet, *Ultrasonics* **42**, 189 (2004).</footnote>

<footnote n="8">⁸ W.J. Xu, F. Jenot, and M. Ourak, *Review of Quantitative Nondestructive Evaluation*. **24**, 156 (2005).</footnote>

<footnote n="9">⁹ DV. Lancellotti and R. Orta, *J. Acoust. Soc. Am.* **104**, 2638 (1998).</footnote>

<footnote n="10">¹⁰ V. Pagneux and A. Maurel, *J. Acoust. Soc. Am.* **110**, 1307 (2001).</footnote>

<footnote n="11">¹¹ O. Bou Matar, N. Gasmi, H. Zhou, M. Goueygou, and A. Talbi, *J. Acoust. Soc. Am.* **133**, 1415 (2013).</footnote>

<footnote n="12">¹² L. Wang and S.I. Rokhlin, *J. Mech. Phys. Solids* **52**, 2473 (2004).</footnote>

<footnote n="13">¹³ C. Baron and S. Naili, *J. Acoust. Soc. Am.* **127**(3), 1307 (2010).</footnote>

<footnote n="14">¹⁴ C. Baron, *Ultrasonics* **51**, 123 (2011).</footnote>

<footnote n="15">¹⁵ C. Baron, *Ultrasound Med. Biol.* **38**(6), 972 (2012).</footnote>

<footnote n="16">¹⁶ A.N. Stroh, *J. Math. Phys.* **41**, 77 (1962).</footnote>

<footnote n="17">¹⁷ J.E. Gubernatis and A.A. Maradudin, *Wave motion* **9**, 111 (1987).</footnote>
to automatically incorporate the mechanical and/or electrical boundary conditions into the equations of motion has been given and explained. Shown are the results for several applications for propagative and stationary waves. Extensions and possible further developments have been addressed.

Position-dependent physical constants are an interesting device which allows to easily incorporate boundary conditions and mechanical and electrical sources into field equations. In addition, they allow to describe exactly heterogeneous and/or inhomogeneous structures.

It is worth mentioning that, for guided waves in a layered plate supported or not by a substrate, the polynomial approach in Ref. 11 and the mapped orthogonal functions method present many similarities. Both use a polynomial expansion, result through an orthogonalization process in a linear eigenvalue problem, can directly determine complex wave numbers for a given frequency and can deal with complicated acoustic wave-based devices, for instance piezoelectric-piezomagnetic composite devices.<fnref n="11" />,⁴⁹ The difference resides in the treatment of the continuity and boundary conditions. For the mapped orthogonal functions method, the continuity of the generalized displacement vector (mechanical displacement components, electric and magnetic potentials) is ensured by a wise choice of the associated polynomial expansion⁴⁷ and the continuity and boundary conditions of the generalized stress vector (normal stress components, normal electric displacement and magnetic induction) are ensured thanks to the δ-functions given by the position-dependent physical constants. All these continuity and boundary conditions are incorporated into equations of motion and the orthogonalization process straightforwardly results in a system with as many equations as unknowns. For (ω real, k complex), the dimensionality of the associated matrix eigenvalue equation for a piezoelectric structure with N<sub>c</sub> homogeneous layers is [2 × N<sub>c</sub> × {3 × (M + 1)}] × [2 × N<sub>c</sub> × {3 × (M + 1)}] where (M + 1) denotes the number of polynomials kept in the expansions of the four components of the generalized displacement vector. For the polynomial approach in Ref. 11 a first part of the system of equations which govern the wave motion is given by the equations imposing the continuity and boundary conditions without orthogonalization process and the second part, through an orthogonalization process, by the equations governing the generalized displacement vector. For (ω real, k complex), the dimensionality of the associated matrix eigenvalue equation for the same piezoelectric structure with N<sub>c</sub> homogeneous layers and the same polynomial expansions is [2 × 4 × N<sub>c</sub> × (M + 1)] × [2 × 4 × N<sub>c</sub> × (M + 1)]. For the eigenvalue problem, the dimensionality is slightly at the advantage of the mapped orthogonal functions method by a factor of nine-sixteenth. The advantage is even more important for piezoelectric-piezomagnetic structures, factor of nine-twenty fifth always at the advantage of the mapped orthogonal functions method. This is counterbalanced, for the approach in Ref. 11, by a lighter preparatory mathematical work when working with more than three distinct sets of polynomials. Probably the polynomial approach in Ref. 11 could be easily extended (i) to deal with through-thickness FGM structures and (ii) to give access, for structures with stationary waves, to modal and forced responses of acoustic wave-based multilayered devices (homogeneous layers supported or not by a substrate).

<footnote n="1">¹ W.T. Thomson, *J. Appl. Phys.* **21**, 89 (1950).</footnote>

<footnote n="2">² L Knopoff, Bull. Seismol. Soc. Am. **54**, 431 (1964).</footnote>

<footnote n="3">³ M.J.S. Lowe, *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **42**(4), 525 (1995).</footnote>

<footnote n="4">⁴ M. Castaings and B. Hosten, in *Ultrasonics International Conference* (Butterworth-Heinemann, Ltd, London, UK, 1993), pp. 431–434.</footnote>

<footnote n="5">⁵ B. Hosten and M. Castaings, *J. Acoust. Soc. Am.* **94**, 1488 (1993).</footnote>

<footnote n="6">⁶ L. Wang and S.I. Rokhlin, *Ultrasonics* **39**, 413 (2001).</footnote>

<footnote n="7">⁷ B. Collet, *Ultrasonics* **42**, 189 (2004).</footnote>

<footnote n="8">⁸ W.J. Xu, F. Jenot, and M. Ourak, *Review of Quantitative Nondestructive Evaluation*. **24**, 156 (2005).</footnote>

<footnote n="9">⁹ DV. Lancellotti and R. Orta, *J. Acoust. Soc. Am.* **104**, 2638 (1998).</footnote>

<footnote n="10">¹⁰ V. Pagneux and A. Maurel, *J. Acoust. Soc. Am.* **110**, 1307 (2001).</footnote>

<footnote n="11">¹¹ O. Bou Matar, N. Gasmi, H. Zhou, M. Goueygou, and A. Talbi, *J. Acoust. Soc. Am.* **133**, 1415 (2013).</footnote>

<footnote n="12">¹² L. Wang and S.I. Rokhlin, *J. Mech. Phys. Solids* **52**, 2473 (2004).</footnote>

<footnote n="13">¹³ C. Baron and S. Naili, *J. Acoust. Soc. Am.* **127**(3), 1307 (2010).</footnote>

<footnote n="14">¹⁴ C. Baron, *Ultrasonics* **51**, 123 (2011).</footnote>

<footnote n="15">¹⁵ C. Baron, *Ultrasound Med. Biol.* **38**(6), 972 (2012).</footnote>

<footnote n="16">¹⁶ A.N. Stroh, *J. Math. Phys.* **41**, 77 (1962).</footnote>

<footnote n="17">¹⁷ J.E. Gubernatis and A.A. Maradudin, *Wave motion* **9**, 111 (1987).</footnote>
18 J.E. Lefebvre, V. Zhang, J. Gazalet, T. Gryba, and V. Sadaune, *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **48**(5), 1332 (2001).

19 L. Elmaimouni, J.E. Lefebvre, V. Zhang, and T. Gryba, *NDT & E Int.* **38**(5), 344 (2005).

20 J.G. Yu, F.E. Ratolojanahary, and J.E. Lefebvre, *Compos. Struct.* **93**(11), 2671 (2011).

21 J.G. Yu, J.E. Lefebvre, C. Zhang, and F.E. Ratolojanahary, *Meccanica* **50**(1), 109 (2015).

22 P.E. Lagasse, *J. Acoust. Soc. Am.* **53**, 1116 (1973).

23 B. Aalami, *J. Appl. Mech.* **40**, 1067 (1973).

24 I. Bartoli, A. Marzani, F.L. di Scalea, and E. Viola, *J. Sound Vib.* **295**, 685 (2006).

25 G. Tobolka, *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **SU-26**, 426 (1979).

26 S.V. Biryukov, Y.V. Gulyaev, V.V. Krylov, and V.P. Plessky, *Surface acoustic waves in inhomogeneous media* (Springer, Berlin, 1995).

27 E.L. Adler, *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **37**(6), 485 (1990).

28 W.R. Smith, H.M. Gerals, J.H. Collins, T.M. Reeder, and H.J. Shaw, *IEEE Trans. Microwave Theory Tech.* **17**, 856 (1969).

29 K.M. Lakin, J. Belsick, J.F. McDonald, and K.T. McCarron, in *Proceedings of the IEEE Ultrasonics Symposium, Atlanta, Ga, USA, October 2001*. pp. 827–831.

30 S. Datta and B.J. Hunsinger, *J. Appl. Phys.* **50**(5), 3370 (1979).

31 A. Raherison, A. J.E. Lefebvre, F.E. Ratolojanahary, L. Elmaimouni, and T. Gryba, *J. Appl. Phys.* **108**, 104904 (2010).

32 G. Endoh, K-Y Hashimoto, and M. Yamaguchi, *Jpn. Appl. Phys.* **34**(5B), 2638 (1995).

33 T. Makkonen, A. Holappa, J. Ellä, and M.M. Salomaa, *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **48**(5), 1241 (2001).

34 S. Datta and Bill J. Hunsinger, *J. Appl. Phys.* **49**(2), 475 (1978).

35 B. A. Auld, *Acoustic fields and waves in solids*, 2nd ed. (Krieger publishing company, Malabar Florida, 1990), Vol. I, p. 43.

36 A.H. Fahmy and E.L. Adler, *Appl. Phys. Lett.* **22**(10), 495 (1973).

37 A.A. Maradudin, R.F. Wallis, D.L. Mills, and R.L. Ballard, *Phys. Rev. B* **6**(4), 1106 (1972).

38 Y. Kim and W.D. Hunt, *J. Appl. Phys.* **68**(10), 4993 (1990).

39 A.A. Maradudin, *Jpn. J. Appl. Phys. Suppl. 2*, 871 (1974).

40 B. A. Auld, *Acoustic fields and waves in solids*, 1st ed. (Wiley-Interscience, New York, 1973), Vol. II, p. 355.

41 J.E. Lefebvre, V. Zhang, J. Gazalet, and T. Gryba, *J. Appl. Phys.* **83**(1), 28 (1998).

42 L. Elmaimouni, J.E. Lefebvre, F.E. Ratolojanahary, A. Raherison, T. Gryba, and J. Carlier, *Wave motion* **48**(1), 93 (2011).

43 J.E. Lefebvre, V. Zhang, J. Gazalet, and T. Gryba, *J. Appl. Phys.* **85**, 3419 (2009).

44 A. Raherison, J.E. Lefebvre, F.E. Ratolojanahary, L. Elmaimouni, and T. Gryba, *J. Appl. Phys.* **108**, 104904 (2010).

45 R.L. Jungerman, P. Bennett, A.R. Selfridge, B.T. Khuri-Yacub, and G.S. Kino, *J. Acoust. Soc. Am.* **76**, 516 (1984).

46 J. Yang, *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **54**(6), 1154 (2007).

47 J.G. Yu, J.E. Lefebvre, and Y.Q. Guo, *Compos. Pt. B-Eng.* **51**, 260 (2013).

48 P.M. Rabotovao, F.E. Ratolojanahary, J.E. Lefebvre, A. Raherison, L. Elmaimouni, T. Gryba, and J.G. Yu, *J. Appl. Phys.* **114**, 124502 (2013).

49 J.G. Yu, J.E. Lefebvre, and C. Zhang, *Compos. Struct.* **116**, 336 (2014).
18 J.E. Lefebvre, V. Zhang, J. Gazalet, T. Gryba, and V. Sadaune, *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **48**(5), 1332 (2001).

19 L. Elmaimouni, J.E. Lefebvre, V. Zhang, and T. Gryba, *NDT & E Int.* **38**(5), 344 (2005).

20 J.G. Yu, F.E. Ratolojanahary, and J.E. Lefebvre, *Compos. Struct.* **93**(11), 2671 (2011).

21 J.G. Yu, J.E. Lefebvre, C. Zhang, and F.E. Ratolojanahary, *Meccanica* **50**(1), 109 (2015).

22 P.E. Lagasse, *J. Acoust. Soc. Am.* **53**, 1116 (1973).

23 B. Aalami, *J. Appl. Mech.* **40**, 1067 (1973).

24 I. Bartoli, A. Marzani, F.L. di Scalea, and E. Viola, *J. Sound Vib.* **295**, 685 (2006).

25 G. Tobolka, *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **SU-26**, 426 (1979).

26 S.V. Biryukov, Y.V. Gulyaev, V.V. Krylov, and V.P. Plessky, *Surface acoustic waves in inhomogeneous media* (Springer, Berlin, 1995).

27 E.L. Adler, *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **37**(6), 485 (1990).

28 W.R. Smith, H.M. Gerals, J.H. Collins, T.M. Reeder, and H.J. Shaw, *IEEE Trans. Microwave Theory Tech.* **17**, 856 (1969).

29 K.M. Lakin, J. Belsick, J.F. McDonald, and K.T. McCarron, in *Proceedings of the IEEE Ultrasonics Symposium, Atlanta, Ga, USA, October 2001*. pp. 827–831.

30 S. Datta and B.J. Hunsinger, *J. Appl. Phys.* **50**(5), 3370 (1979).

31 A. Raherison, A. J.E. Lefebvre, F.E. Ratolojanahary, L. Elmaimouni, and T. Gryba, *J. Appl. Phys.* **108**, 104904 (2010).

32 G. Endoh, K-Y Hashimoto, and M. Yamaguchi, *Jpn. Appl. Phys.* **34**(5B), 2638 (1995).

33 T. Makkonen, A. Holappa, J. Ellä, and M.M. Salomaa, *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **48**(5), 1241 (2001).

34 S. Datta and Bill J. Hunsinger, *J. Appl. Phys.* **49**(2), 475 (1978).

35 B. A. Auld, *Acoustic fields and waves in solids*, 2nd ed. (Krieger publishing company, Malabar Florida, 1990), Vol. I, p. 43.

36 A.H. Fahmy and E.L. Adler, *Appl. Phys. Lett.* **22**(10), 495 (1973).

37 A.A. Maradudin, R.F. Wallis, D.L. Mills, and R.L. Ballard, *Phys. Rev. B* **6**(4), 1106 (1972).

38 Y. Kim and W.D. Hunt, *J. Appl. Phys.* **68**(10), 4993 (1990).

39 A.A. Maradudin, *Jpn. J. Appl. Phys. Suppl. 2*, 871 (1974).

40 B. A. Auld, *Acoustic fields and waves in solids*, 1st ed. (Wiley-Interscience, New York, 1973), Vol. II, p. 355.

41 J.E. Lefebvre, V. Zhang, J. Gazalet, and T. Gryba, *J. Appl. Phys.* **83**(1), 28 (1998).

42 L. Elmaimouni, J.E. Lefebvre, F.E. Ratolojanahary, A. Raherison, T. Gryba, and J. Carlier, *Wave motion* **48**(1), 93 (2011).

43 J.E. Lefebvre, V. Zhang, J. Gazalet, and T. Gryba, *J. Appl. Phys.* **85**, 3419 (2009).

44 A. Raherison, J.E. Lefebvre, F.E. Ratolojanahary, L. Elmaimouni, and T. Gryba, *J. Appl. Phys.* **108**, 104904 (2010).

45 R.L. Jungerman, P. Bennett, A.R. Selfridge, B.T. Khuri-Yacub, and G.S. Kino, *J. Acoust. Soc. Am.* **76**, 516 (1984).

46 J. Yang, *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **54**(6), 1154 (2007).

47 J.G. Yu, J.E. Lefebvre, and Y.Q. Guo, *Compos. Pt. B-Eng.* **51**, 260 (2013).

48 P.M Rabotovao, F.E. Ratolojanahary, J.E. Lefebvre, A. Raherison, L. Elmaimouni, T. Gryba, and J.G. Yu, *J. Appl. Phys.* **114**, 124502 (2013).

49 J.G. Yu, J.E. Lefebvre, and C. Zhang, *Compos. Struct.* **116**, 336 (2014).