Finite Element Model for SAW

In contrast to the well known Rayleigh waves which consists of partial lon-gitudinal waves and shear waves, the LSAW mainly propagates in the shear direction on the sagittal plane and is trapped at substrate surface and satis-fies the stress free boundary condition on the surface. These properties allow one to reduce the general mode analysis in 3D to a 2D problem as shown in Figure 1(b) [28]. Furthermore, the boundary conditions for displacement can naturally be set to be rigid on the bottom boundary and stress-free on the top surface, and the boundary conditions for the electric potential can be set to be short-circuited for the electrodes on the top boundary and open-circuited else-where [10]. As proved in Auld’s book [4], these boundary conditions guarantee the mode orthogonality and further ensure the mode excitation is determined

(a) A standard configuration of SAW res-onators

(b) A 2D model of a LSAW resonator on the sagittal plane

Figure 1:

by the applied traction force and potential on the free surface. Therefore, on the sagittal plane, the usual 2D mode analysis can be applied to analyze the LSAW on the resonators with IDTs. In the following, we only consider the LSAW resonator on a 2D plane (the sagittal plane associated with crystal cuts).

To model the wave propagation in an infinite domain with periodically ar-ranged electrodes, thanks to the Floquet-Bloch Theorem, one can reduce the problem to a single cell domain with one IDT by assuming the wave ψ is quasi-periodic of the form ψ(x1, x2) = ψp(x1, x2)e^(α+iβ)x1 where x1is the wave prop-agation direction, p is the length of the unit cell (i.e. the periodic interval), α and β are the attenuation and phase shifts along the wave propagation direc-tion, respectively, and ψp satisfies ψp(x1+ p, x2) = ψp(x1, x2). Let Ω denote the PZT with a single IDT as shown in Figure 2, and Γland Γr denote the left and right boundary segments of Ω. For the general anisotropic PZT substrates, under the assumption of linear piezoelectric coupling, the elastic and electric fields interact following the general material constitutions below

T = c^ES− e^⊤E,

D = eS + ε^SE, (3)

where vectors T , S, D and E are the mechanical stress, strain, dielectric dis-placement and the electric field, respectively, and the matrices c^E, ε^S and e are the elasticity constant, dielectric constant and piezoelectric constant matrices measured at constant electric and constant strain fields at constant temperature.

For various crystal cut of the PZT, the material constant matrices c^E, ε^S and e depend on the Euler angle θ of the cut. By applying the Bond strain trans-formation matrix N_θ [5] and the usual coordinates transformation matrix M_θ to the strain field and electric field, respectively, the material constant matrices for the cut angle θ can be obtained by

c^E := [Nθ]^⊤c^E₀[Nθ], e := [Mθ]^⊤e0[Nθ], and, ε^S := [Mθ]^⊤ε^S₀[Mθ],

Figure 2: A 2D single cell domain of a LSAW resonator and boundary conditions

here c^E₀, e0, and ε^S₀ denote the material constant matrices of the crystal cut at Euler angle θ = [0^o, 0^o, 0^o].

By applying the virtual work principle to the equation (3), the equilibrium state under the external body force f , the electrical field g and the above men-tioned boundary conditions of the LSAW resonator, we have

∫ displacement vector, ϕ is the electric potential that satisfies ∇ϕ = E, S = [^∂u_∂x,^∂v_∂y,^∂w_∂z,^∂v_∂z+^∂w_∂y,^∂w_∂x+^∂u_∂z,^∂u_∂y+^∂v_∂x]^⊤, and δq, δϕ and δS are the corresponding virtual displacement, potential and strain vectors, respectively. The equation can then be discretized by finite element method [1, 10]. Following the usual free mode analysis, we consider f = 0, g = 0 and a time harmonic quasi periodic solution vector ψ_ω(x, t) = ψ(x)e^iωt. The spatial function ψ(x) = [q(x), ϕ(x)] sat-isfies the boundary conditions shown in Figure 2 in which the periodic boundary

conditions, proposed by Buchner [6], (4), the equation can be rewritten in the following matrix form:

[ K^qq− ω²M^qq K^qϕ Mechanical damping effects can also be considered by using the Rayleigh

damp-ing assumption in which the matrix K^qq − ω²M^qq in (8) are modified into K^qq+ iω(κ1K^qq+ κ2M^qq)− ω²M^qq. Here κ1and κ2are coefficients associated with the viscous damping and mass damping, respectively.

To obtain the palindromic quadratic eigenvalue problem associated with the propagation parameter γ, following Hofer’s approach [14], the nodal unknowns are splitted into the inner nodes ψ_i = [q_i, ϕ_i], the left boundary nodes ψ_l = [q_l, ϕ_l] and the right boundary nodes ψ_r= [q_r, ϕ_r]. The matrix equation (8) can be recasted into the following form:



here Rland Rrare vectors obtained from the discretization of the terms Fl+ Ql

and Fr+ Qr, respectively. From the periodic boundary conditions (5), (6) and (7), (9) becomes

Furthermore, by multiplying the matrix [ Ii 0 0

0 γIl Il

]

to (10), the GEP associated with the propagation parameter γ is obtained:

([ K_ii K_il

Since the viscosity is small for most crystalline solids, the attenuation factor α is close to zero. As a result, the propagation factors λ near the unit circle, denoted byU, is desired. Moreover, for the frequency ω in the stopping band, the frequency shift parameter β shall be close to π when the periodic interval p (i.e. the domain width here) equals to half of the incident wave length λ₀. Therefore, for the eigenvalue problem (11), we are interesting in finding the eigenvalues λ close toU, especially for those are near −1 on the complex plane.

Notice that the nonzero eigenvalues of (11) appear in the reciprocal pairs (λ, 1/λ). The reciprocal relation is very sensitive to numerical errors when they are close toU. On the other hand, it is well known that the solution of a general elliptic problem have singularities around corners [12] and, in addition, the solu-tion may become less regular near the interface between the electrode and PZT substrate. It is inevitable that the error from discretization may be amplified in computing the reciprocal pairs. Therefore, it is important to minimize the accuracy deterioration due to singularities and lower regularity in finite element solutions. One can resolve the singularity by constructing the singular elements in which the mesh points are clustered to the singular source according to the order of the singularity [3]. In our calculation, we simply employee the locally refined meshes. An additional benefit from using the locally refined meshes is that we can discretize equation (8) using linear elements instead of using high order finite element discretization [6]. However, drawbacks include (i) the ma-trices from the discretization of (8) is large and sparse and (ii) the sparse pattern of the matrices is unstructured. These make the efficient computation of eigen-values for large GEP in (11) a challenge. Moreover, for pizoelectric crystals, the elastic constant matrix c^E is 10²⁰ greater than the electric constant matrix ε^S. To compute the eigenvalues and eigenvectors, proper scaling between the mechanic field and the electrical field is required. Hence, the eigen solutions ob-tained from the scaled problem must be accurate enough in order to disregard the round off error in the re-scaling process. Therefore, for solving the large sparse eigenvalue problem (11), an efficient algorithm, not only preserves the reciprocal eigen-structure but is also accurate enough to prevent error amplifi-cation from rounding and discretization, is desired. In the next section, we shall introduce an efficient structure-preserving algorithm that ensure the accuracy of the eigen-curves λ(ω) and the associated eigenvectors of (11).