As stated, this section is devoted to the resolution of the mechanical equations that govern the PMUT plate, in air and, when the plate is coupled with a fluid medium. Mechanical equations are first developed, described, and then solved by means of a FD discretization scheme. Second, a first model validation will be presented by comparing the simulated resonance frequencies with those provided by a commercial finite element code (COMSOL Multiphysics [
54,
55]). Then, the fluid-PMUT coupling equations will be introduced in numerical form, after the FD discretization step. The last subsection is the second validation step of the model, where simulated results are compared with analytical models from the literature.
Before providing the mechanical plate description, the typical geometry of a PMUT cell with the corresponding coordinate system axis is established. Our PMUT structure plate is made up of three layers (
Figure 1) and is axisymmetric. The layers include a structural layer, usually silicon, a piezoelectric layer, and a partially covering top electrode. Here to be consistent with the experimental devices presented in the last section, there is no metallic bottom electrode; because the silicon material is doped, it also plays the role of an electrode. The plate is considered to be mechanically clamped around its outer perimeter. In accordance with the PMUT circular geometry, the plate analysis will be carried out in cylindrical coordinates (
,
,
). The terms
and
are the plate radius and the top electrode radius, respectively. The position of the top of the
k-indexed layer along the z-axis is referenced with the
coordinate, where
, and (
-) corresponds to the thickness of the
k-layer. Note that to clearly identify the piezoelectric layer from the others, its upper coordinate is labeled
. The dashed line with coordinate
is the reference line, the so-called neutral line, whose expression and value will be given later.
2.1. Mechanical Behavioral Equations and Resolution With Finite Difference Discretization
The mechanical behavior of the PMUT plate is based on the resolution of the classical Kirchhoff-Love thin plates equation [
56,
57,
58] expressed under cylindrical coordinates [
46,
48,
59,
60]. This theory relies on two mains assumptions:
- 1)
The mechanical plate vibration is limited to the displacements (,) and (,) along the r-axis (axial displacement) and -axis (transverse displacement), respectively.
- 2)
The through-the-thickness stresses and strains are negligible.
Total displacements
and
are shown in
Figure 2. Using the physical neutral line concept and classical plate theory, the displacements take the following forms [
47,
58]
where
is the axial displacement of the reference line
,
i.e, the neutral line, which will not be considered negligible. In many papers [
35,
47,
61,
62], authors deliberately neglect this axial displacement, which enables to strongly simplify the analytical calculations. However, in papers from Sammoura’s group [
49,
63], it was shown that the axial displacement of the plate plays a role in the piezoelectric equation of the PMUT and cannot be suppressed. Therefore, at this step of the theoretical developments, we chose to keep this axial displacement in the set equations and later discuss its influence on the PMUT behavior. Due to axisymmetric deformation, only the strains
and
must be considered, and the shear strain
is null. Using the set of relations (1) and (2),
and
are expressed as follows:
To establish the dynamic equilibrium equation, we classically introduce the resultant stresses that apply to the plane perpendicular to the neutral line [
56,
57,
58]:
The corresponding resultant bending moments are defined as [
56,
57,
58]:
where
and
are stress field vector components according to the
and
directions, respectively.
To introduce the expression of the displacements (1) in the resultant stresses and resultant moments, one must use the classical strain-stress relation (Hooke’s law) and then compute the integral along the
-axis for each term of
Trr and
Tθθ. These calculation steps have already been widely described in the literature [
58] and will not be detailed here. However, since axial displacement is not suppressed from the equations, we need to introduce the following terms:
- 1)
The membrane strain matrix:
- 2)
The flexural strain matrix, known as the curvature matrix:
The resultant moments expression becomes:
where the
matrix is written as:
The resultant moments expression can be reorganized thanks to the following decomposition:
where
and
are the resultant moments associated with the
and
terms, respectively. If we write equation (9) with matrix operators A, B and C, developed by Zhang
et al. [
47] and by J.N. Reddy’s book [
58], one has:
where
,
,
,
,
and
.
Similar expressions can be obtained for Nrr and Nθθ, and the calculations and relations are given in the Appendix section.
The equations of motion of the plate for circular vibrations, without an external source, are given by J.N. Reddy [
58]; for the resultant moments, we have :
For the resultant stresses, if there are no external forces acting in the volume of the plate, we have:
where
is the plate weight per unit area. Note that relation (11) assumes that the rotary inertia terms are neglected, as explained by J.N. Reddy [
58].
is a source term that corresponds to an external distributed pressure that applied to the plate.
Expressions of
and
are introduced in relation (11)(10) to obtain:
The first three terms on the left-hand side of relation (13) depend on
, and the three others depend on
. If the axial displacement of the neutral line
was neglected, relation (13) would be a function of only the transverse displacement
. A technique to separate transverse and axial displacement components is to modify the expression of these first three terms as follows:
If the neutral line coordinate
zn is chosen such that:
then the terms with
are removed from the relation (13).
The neutral line is the position to uncouple the axial displacement to the transverse displacement and the exact expression of
is then deduced:
We find the expression used in the literature [
49,
62], but it has never been demonstrated or even verified in the case of acoustic MEMS. Only Sammoura proposed a demonstration [
49], but the initial hypothesis of his model considered the displacement
to be negligible. The same result is obtained in the appendix section for the second equation of motion.
The final equation of motion becomes:
where
and
is the Laplacian operator in cylindrical coordinates.
Equations (18) and (19) are the equivalent flexural rigidity and the equivalent Poisson’s coefficient of the plate, respectively.
The solution of equation (17) can be obtained with analytical developments. In previous reports, it is common to decompose the solution into a basis of eigenmodes and to determine the coefficients of this decomposition to meet the boundary conditions. We did not choose this approach because the function basis used for air coupling does not properly account for when the PMUT plate is coupled with a fluid medium (see
Section 2.3).
To solve this equation the FD discretization scheme presented in [
57] was used. Our group previously used this scheme for circular-shaped CMUTs [
53]. To convert equation (17) from its analytical to its numerical form, we first need to set up discretization nodes (see
Figure 3) along the plate: from the center at
r = 0 to the edge of the plate at
r = a. The pitch between two nodes
Δr is fixed, and the number of nodes is labeled N. At each discretization node
i, relation (17) is combined with the discretization scheme given by Timoshenko [
57] to obtain a relation between the displacement
of the
i-node with its neighboring nodes. Particular attention is given to the first node and the last node to take into account the boundary conditions. In our case the following boundary conditions were used:
- (a)
, symmetrical boundary conditions,
- (b)
, clamped boundary conditions.
The set of relations obtained is then gathered to express the relation (17) in matrix form:
is the vector made with all degrees of freedom
; in air there is no external pressure, so [
q] is null.
is the stiffness matrix that results from the discretization of the differential operator of relation (17), and
is the mass matrix.
is a pure diagonal matrix that contains the value of
at each discretization node. The stiffness, mass and neutral line associated with the nodes before and after the boundary are thus perfectly defined.
2.2. Comparison with the Finite Element Model (FEM)
The first validation step was to compare the simulation results of the FD model with those provided by a commercial finite element model. COMSOL Multiphysics was used. There are two key points to check. The validation of the neutral line coordinate
computation in the case where the plate thickness is not homogeneous is the first concern. In other words, is the model able to predict the plate response when the neutral line position is discontinuous? To overcome this difficulty, some authors have proposed separating the plate into two zones and combining the solutions by applying continuity conditions at the interface between the two zones [
61,
64,
65,
66]. However, even if analytical solutions are developed, the model becomes quite "heavy" as it introduces many additional variables.
The anisotropy of the structural layer is the second concern of this validation step. All equations were developed in the case of axisymmetric vibration, assuming that the materials are isotropic. In practice, the structural material is often silicon material, which is anisotropic. As silicon is an anisotropic material [
67], the Young’s modulus may vary from 130 GPa to 188 GPa according to the axis [
] or [
] in crystal lattices. In a previous study [
68], it was shown that a circular isotropic silicon plate with a Young’s modulus of 143 GPa and a Poisson’s ratio of 0.278 gives the same deformation and resonance frequency as a circular anisotropic silicon plate, with less than 0.3% error. To conduct simulations, it was chosen to model silicon material as an isotropic material and to keep these elastic parameter values.
To compare the FD model with COMSOL Multiphysics, we defined a PMUT cell comparable to the literature structures based on Aluminum Nitride (AlN), since the devices used experimentally are also AlN-based (see
Section 4). The PMUT cell is made of three layers of only doped silicon, which is the structural layer and the bottom electrode, an AlN layer, and an Aluminum (Al) layer for the top electrode. The layer thicknesses and PMUT diameter were fixed to obtain a resonance frequency in air and in water in the MHz range. The properties of AlN used were extracted from the literature [
69] and from the experimental results presented in the last section of the paper,
Section 4. All material properties and dimensions are reported in
Table 1, but only the mechanical properties were used for this comparison. For the FD model 50 nodes of discretization were used along the
r-axis. For the FEM model, no asymmetrical boundary conditions were applied to the silicon layer, and the real elastic properties (anisotropic) [
67] were used. The AlN and aluminum layers are not difficult to simulate since AlN has hexagonal crystal symmetry and aluminum is isotropic. The mesh of the studied structure was carefully chosen to meet the λ/4 criterion in the whole frequency range so 125 hexahedron-shaped elements in total were chosen.
Figure 4 shows the eigen-frequency of the first PMUT plate resonance according to the top electrode radius, which varies from 0 to 50 µm (100% surface metallization rate). The FDM predictions perfectly match the FEM simulations. The maximum relative error is 1.13% when the normalized electrode radius is close to 50%. When the electrode radius increases, a decrease in the resonance frequency is observed due to a mass effect; then, over a normalized electrode radius of 80%, the resonance frequency increases. This is ascribed to an increase in the global plate stiffness.
This first step clearly justifies the choice of modeling the silicon layer as an isotropic material and the numerical solving of the plate equation with discontinuity of the neutral line position.
2.3. PMUT/fluid Coupling: Implementation of a Boundary Element Matrix
In this third subsection, the model is modified to introduce the coupling with a fluid, which is modeled through an external distributed pressure
. From the plate equation (20), one term must be added:
is the radiated pressure vector, which depends on the displacement vector
. The two vectors
and
are linked by the radiation boundary conditions at the PMUT/fluid interface. This relation can be expressed mathematically using a boundary matrix
:
where each term of the
matrix,
corresponds to the mutual acoustic impedance between the
i-indexed node and the
j-indexed node. In other words, this term represents the force exerted by node
i when it radiates with displacement
onto node j. The radiating surface associated with node
i is a ring centered at the
position with Δr width and the surface associated with node
j is a ring centered at the
with the same width. If the discretization pitch is thin enough, it is safe to assume each ring vibrates like a piston. To compute the
matrix terms, one has to use the Huygens-Rayleigh integral [
63]:
where
and
are the surfaces of
i-indexed ring and
j-indexed ring respectively.
is Green’s function of the fluid for rigid baffle conditions:
is the acoustic wavenumber in the fluid.
and
are acoustic properties (density and velocity respectively) of the fluid. The numerical computation of the first integral (integration over the emitting ring of surface
) in relation (23) can be extremely time-consuming because the convergence is typically very slow and strong numerical instabilities may be present. Seybert
et al [
70] thus proposed splitting the integral into two parts to stabilize computations and increase the convergence time. We have used this decomposition. The validation of the
matrix terms computation was performed assuming the PMUT vibrates as a piston, i.e, each term of the
vector is equal to 1. The radiation impedance of the circular piston computed numerically with the FD model was compared with the analytical expressions [
71]. No difference was observed.
2.4. Comparison with Literature
For this second validation step, the analysis focuses on PMUT/fluid coupling, where the simulations obtained from the FD models are compared with those obtained with analytical models from the literature. The present goal is to discuss the limits of each approach and identify their respective advantages and disadvantages. Many analytical models have been published, but the two oldest [
12,
46] have been used as the basis for assumptions applied in the more recent works [
33,
59,
61,
72]. We have therefore chosen to focus our comparison on the model from K. Smyth
et al [
46] and Perçin
et al [
4,
12].
Perçin [
12] developed a lumped-element circuit to model PMUT/fluid coupling. From a side, the PMUT plate mechanical impedance is computed when the device is vibrating in air. A uniform external pressure source is applied over the whole plate, and the plate equation (equation (11)) is solved to obtain the average plate velocity and thus the plate mechanical impedance. Analytical resolution of the plate equation based on modal decomposition of the solution is used. From another side, the plate radiation impedance of the PMUT is computed assuming that it vibrates like a piston with rigid baffle boundary conditions. The published radiation impedance expression of a piston in fluid [
71] is directly used. K. Smyth
et al [
46] developed the same approach to compute the plate mechanical impedance. The two models differ for the radiation impedance determination. K. Smyth
et al [
46] took into account the nonuniform velocity distribution along the plate to compute the radiated pressure. The velocity distribution along the plate radius is the same as that obtained from the modal decomposition in air. This approximation is clearly better than considering that PMUT works like a piston, particularly at high frequencies where the velocity distribution impacts the PMUT directivity. Pala
et al [
72] used the same method as Smyth except that a decomposition of quadratic-type functions was used instead of eigenfunctions of the plate equation.
To compare the three models, the device presented in
Section 2.2 with a top electrode radius of 35 µm, so a metallization ratio radius of 70% was simulated. For the FD model, a uniform acoustic pressure of 1 Pa was applied (
i.e, the
vector in equation (21)), and the spatial average plate velocity was computed from the computed displacement field vector
. For the analytical models, we implemented the analytical solutions provided by Smith [
46] and Perçin [
12] and computed the impedance terms. The average velocity is straightforward to determine from the lumped-element equivalent electrical circuit.
Figure 5 shows the spatial average particle velocity of the PMUT according to the frequency for the three models in water, and the simulation of the device working for the FD model in air is added. The resonance frequency in air is close to 2.5 MHz, and it drops to a value close to 1 MHz in water. The three models yield very similar resonance frequency values in water, thereby validating the numerical implementation of the PMUT/fluid coupling in the FD model. However, the FD model differs from the two analytical models at higher frequencies. In air, at 9 MHz, a cutoff frequency appears that corresponds to the situation where the spatial average velocity of the plate is null. In water, with the FD model, we see that this cutoff frequency drops to 5 MHz, while the two analytical models provide 9 MHz, as in air. This cutoff frequency is a key characteristic of the PMUT cell since it also defines the high cutoff frequency of a PMUT array. This frequency remains the same whether the PMUT operates alone or as part of an array. This is clearly one of the limitations of the analytical models presented to date since the mechanical impedance of the PMUT plate is calculated when the device is vibrating in air. However, there are no simple analytical solutions to the plate equation when coupled to a fluid. It is therefore essential for the determination of equivalent mechanical and radiation impedances that the mechanical deformation of the plate used considers the coupling with water. This is the advantage (see
Section 3.2) of having developed a numerical approach. However, it should be remembered from previous works, K. Smyth [
46] and Pala [
72], that when calculating radiation impedances, it is preferable to take plate deformation into account, particularly for the determination of the elementary directivity acoustic field but also the acoustic mutual impedances matrix between PMUTs. Though it is not further in this paper, additional investigations indicated that the type of pressure application and its distribution along the plate impact the cutoff frequency value as well. These additional investigations are ongoing and will be presented in a future.