Simulating the magnetic fields generated by piezoelectric devices using FEM software: Beyond the quasistatic approximation

Piezoelectric materials are anisotropic materials that enable the transfer of energy between acoustic (i.e., elastomechanical) and electromagnetic fields. Hybrid vibrations comprising both electromagnetic and acoustic components can be induced by either driving the material mechanically or by exposing it to an oscillating electric field (most often produced by a pair of electrodes).

One example of electromagnetic-acoustic (EMA) coupling is a surface acoustic wave (SAW) on a piezoelectric substrate. Here, the hybrid wave is localized on the substrate’s surface. Both the acoustic and the electromagnetic parts of such a SAW propagate along the surface with the same frequency and wavenumber. Each part is uniform in the transverse direction. Devices exploiting piezoelectric SAWs encompass delay lines, filters, and resonators; such SAWs also provide a way of coupling phonons to artificial atoms¹ in quantum systems.

Another example of EMA coupling is a piezoelectric antenna. To radiate efficiently, the length of a conventional antenna needs to approach the wavelength of the electromagnetic wave being emitted. Because elastomechanical waves propagate much more slowly (by a factor of around $104$⁠), the size of a piezoelectric antenna operating at the same frequency is vastly smaller. If such an antenna could be gotten to radiate at reasonable efficiency (electromagnetically), it would offer huge potential for the further miniaturization of telecom devices.

Owing to the anisotropy of the medium, the complete simulation of EMA coupling in piezoelectric devices can rarely be done analytically. The so-called quasi-static approximation is often applied to simplify the set of equations to be solved. Under this approximation, the magnetic field is completely neglected, and the electric field strength is represented by the electric potential $ϕ$ as $E_=−∇ϕ$⁠. In many piezoelectric applications, this approximation works well because the magnetic field stores far less energy than what the electric field does.² However, in certain applications, such as antennas³ or high magnetic-Purcell-factor resonators,⁴ knowledge of the magnetic field is crucial: one cannot calculate the antenna’s efficiency (radiation impedance) without it. The quasi-static approximation must be avoided.

Certain known analytical methods can precisely simulate the EMA coupling without invoking the quasi-static approximation. Examples include the superposition-of-partial-waves method,² transverse-resonance method,⁵ and Kirchhoff-plate-theory method.^6,7 These methods are limited to very specific spatial geometries and crystal classes (point group symmetry)—they cannot be extended to simulate arbitrary 3D structures or piezoelectric materials of a different class. Alas, many it not all of the popular finite-element-method (FEM) based simulation platforms, like COMSOL Multiphysics (both its structural mechanics and MEMs modules), implicitly invoke the quasi-static approximation, meaning that they cannot be used to simulate the complete electromagnetic behavior of piezo-mechanical devices. The prime purpose of this paper is to remedy this frustrating inadequacy. Our open method is applicable to any FEM simulation software where weak forms define the whole system of coupled partial differential equations. COMSOL is used here to provide explicit examples. Results for a Rayleigh mode propagating in the X-direction on the Z-cut surface of a $LiNbO3$ crystal are compared against both experimental results and analytical solutions, so as to validate our method’s correctness and accuracy. The results corresponding to the torsional mode on a cylindrical piezoelectric antenna showcase the capability of our method to simulate complicated devices. Definitions of variables, example codes and various details and tips relevant to the building of working simulations are provided in the supplementary material.

The total energy density $U$ at an arbitrary point in space is represented using 8 variables: velocity $v_$⁠, momentum density $p_$⁠, strain $S__$⁠, stress $T__$⁠, electric field strength $E_$⁠, electric displacement $D_$⁠, magnetic field strength $H_$ and magnetic flux density $B_$⁠,² 

where “:” denotes the double dot product between tensors. Note that this form of $U$ restricts our analysis to linear materials.

The above weak form is written as an integration of the negative total energy density throughout space, where one of the two variables in each term is replaced with a test function (denoted by “$∼$”),

In a piezoelectric material, the coupling between acoustic fields and electromagnetic fields is described using the constitutive relations

where $c____E$ is the stiffness tensor at zero or constant electric field strength, $e___$ and $e___T$ are the piezoelectric coupling tensor and its transpose, and $ε__S$ is the permittivity matrix at zero or constant strain. The number of underlines (“$_$”) denotes a tensor’s rank.

For simplicity, we assume that the piezoelectric material is non-magnetic, i.e.,

though extending our analysis to linear magnetic media with non-unity relative permeabilities would be wholly straightforward. We use Faraday’s law to connect the electric field strength to the magnetic flux density

Here, (3a), (3b), (4), and (5) enable (2) to be represented by only three variables $ω$⁠, $u_$⁠, and $E_$⁠. In COMSOL, the last two variables are named Dependent Variables, which are defined when establishing a new model. The components of (2) then serve as the inputs into Weak Form PDE and Weak Contribution nodes in COMSOL’s various different physics interface(s). The independent variable angular, frequency $ω$⁠, is finally solved in Eigenfrequency study. The detailed, explicit representation of our method (defining how each contraction is performed—for example) is provided in the supplementary material.

Like Ref. 8, a penalty term is used to regulate the electromagnetic field. This penalty term is modified from Gauss’s Law. Since the charge density in the device is zero, the penalty term is written as

As mentioned in Sec. II A, the components of electric field strength ${Ex,Ey,Ez}$ are usually chosen as the dependent variables for the electromagnetic field. The subsequent derivation of the magnetic field strength, however, requires division by the angular frequency $ω$ (of the simulated mode), which is a priori unknown. This paradoxical situation can be gotten around by constructing ${Mx,My,Mz}$ as the dependent variables. These variables are defined as

The performance of a piezoelectric structure as a magnetic-field generator can be quantified by its magnetic coupling coefficient, defined as

where $Ekin$⁠, $Eela$⁠, $Eelec$⁠, and $Emag$ denote kinetic energy, elastic energy, electric energy, and magnetic energy, respectively. The superscripts “$int$” and “$ext$” denote the internal field (within the piezoelectric substrate) and external field (in the surrounding environment). Each energy is calculated using an integration of corresponding energy density throughout its area.

Two test cases are illustrated below to verify the accuracy and capability of our method. These cases correspond to two of the most popular applications of saw devices: SAW waveguides and piezoelectric antennas.

We here purely focus on how the SAW propagates—as opposed to how is it generated and detected using interdigital transducers (IDTs). The domain analyzed can be divided into two cuboids: one contains the piezoelectric substrate; the other contains the surrounding environment. These two cuboids are separated by the surface of the substrate. Spatial coordinate axes are chosen such that $x$ is oriented along the direction of propagation, $y$ lies along the surface and is perpendicular to $x$⁠, while $z$ is outward normal to the surface. When no rotation is applied, the material coordinates (or crystal axes used by Auld²) align the spatial coordinates.

Compared to others, a Rayleigh mode generates a relatively strong and moreover analytically calculable magnetic field. The verification of our method thus focuses on its simulation. A schematic diagram of the Rayleigh mode is shown in Fig. 1. Boundary conditions suitable for this mode are applied: A periodic boundary condition is applied on the propagating direction [surfaces ① and ② and their opposites in Fig. 2(a)]. The distance between these two surfaces is set as $100μ$m, and this distance will become the wavelength of the Rayleigh mode. Another periodic boundary condition is applied in the $y$ direction [surfaces ④ and ⑤ and their opposites in Fig. 2(a)]. Since the Rayleigh mode is uniform in this direction, the distance between the two surfaces does not influence the form of the solution. This distance is set as $10μ$m. The fixed constrain is applied in the $z$ direction [surface ③ and its opposite in Fig. 2(a)]. The Rayleigh mode is evanescent in this direction so that $300μ$m, which is three times as long as the wavelength, is set as the distance from the surface to the fixed constraint. The electromagnetic continuity boundary condition is applied on surface ⑥, whereas this boundary is set as a free boundary for the acoustic field. For some metallized surfaces, an electric-wall boundary condition needs to be applied on surface ⑥. In this simulation, $z<0$ [Area $A$ in Fig. 2(a)] is filled with piezoelectric materials and $z>0$ [Area $B$ in Fig. 2(a)] is chosen to be air. $z=0$ is the surface of the substrate. COMSOL’s built-in material properties are used.

The acoustic displacement field, electric field, and magnetic field of the X-propagating Rayleigh mode on Z-cut $LiNbO3$ are illustrated in Figs. 2(b)–2(d), respectively; the fields shown have been normalized. Comparison among experimental results, analytical results, and our FEM method results of SAW velocities is shown in Table I. Velocities of different materials and different modes with different boundaries are included. The relative errors of our results compared with analytical results are defined as

where $v$ denotes our results and $vA$ denotes the analytical results. The analytical results are calculated by ourselves using the method introduced by Tseng.⁹ All results of our method were calculated using 1.7 $μ$m element size. Note that the built-in densities of CdS and ZnO are different from standard values so that the standard values should be input manually. Furthermore, the slowness surface of $128°$ Y-cut $LiNbO3$ is simulated and plotted (Fig. 3). Experimental results from the literature¹⁰ are also included for comparison. We find that our FEM-simulation results lie very close in value to the analytically calculated ones. Both lie reasonably close to what is measured experimentally. These proximities strongly suggest that our method is correct and accurate. The small residual differences stem from experimental measurement errors, the choice of tolerance in the bi-section method used to solve the analytical equation, and the fineness of the mesh used in our FEM method.

The torsional mode also has a significant magnetic field. The material is exactly on the equatorial plane and the z-axis is motionless. Beyond this nodal plane and axis, the material rotates around the z-axis, the amplitude of displacement increases with the distance from the origin [shown in Fig. 5(1)]. In uni-axial materials like PZT-5H, this mode is significant. A cylinder made with other materials, like $LiNbO3$⁠, will show a torsional mode only when the aspect ratio (defined as the ratio of height to the radius $h/r$⁠) is substantial. [There is no hard cut-off but we would suggest $h/r>5$⁠.] If the aspect ratio is small, the frequency of the torsional mode is very high, with many other competing modes close to it (in frequency) (Fig. 4).

Simulations were conducted on a series of PZT-5H cylinders, varying in height from 1 to 5 mm, all with a fixed radius of 10 mm. In each simulation, the cylinder resides at the center of a sphere, representing the surrounding environment. A continuous boundary condition is applied to the cylinder’s surfaces, while, an absorbing boundary condition is applied on the sphere’s outer surface. The acoustic displacement field, electric field, and magnetic field on the central YZ cross section are shown in Fig. 5 as an example. It is found that when the aspect ratio decreases, both the magnetic coupling coefficient and frequency will increase. The relation among those three quantities is shown in Fig. 6. The highest coupling coefficient in this model is $2.6×10−7$ for a cylinder with a small aspect ratio of (⁠$h/r=0.1$⁠). In reality, the torsional mode is unlikely to appear on a disk with a very small aspect ratio so the magnetic coupling coefficient should not be very high in this system.

Verification of the accuracy of our model was performed on fundamental modes. Shaw¹⁴ measured the frequency of axial symmetric normal modes on barium titanate disks with different aspect ratios. The same modes were simulated using our method and the corresponding resonant frequencies and acoustic displacement fields were compared with Shaw’s results. A barium titanate cylinder was used in the simulation, whose height was set as 5 mm and radius changed from 2.5 to 17.5 mm. The properties (elasticity matrix, coupling matrix, and relative permittivity) of barium titanate ceramic were obtained from Dent.¹⁵ The simulated and measured frequencies are plotted in Fig. 7. The simulated acoustic displacement fields at the upper surface of the disk with two different aspect ratios are plotted and compared in Fig. 8. These proximities clearly suggest that our method is correct and accurate. The discrepancies stem from experimental measurement errors, the fineness of the mesh used in our method, and the accuracy of the properties of barium titanate. Note that an inverse aspect ratio $2r/h$ was used in Shaw’s paper. For consistency, we use the same quantity as the independent variable.

This paper fills a fundamental “blind spot” in the electromagnetic simulation of devices involving piezoelectric crystals/substrates: mechanical motion in piezoelectric materials generates magnetic fields, which need to be precisely quantified. A lack of such simulation leads to numerous limitations like one cannot calculate the Poynting vector, which is the radiation energy flow of an antenna. To remedy the problem, we have constructed an efficient, cost/resource-effective simulation method, as described in this paper. This method is generic, based on the so-called “weak forms” and can be implemented in any finite-element-method simulation software. Because of its fundamental scope and the many areas in applied device physics it connects to, we believe our handy method will be useful in the simulation of many fields, such as SAW-based resonators, piezoelectric antennas, and magnetically sensitive quantum devices involving spins.

See the supplementary material for an introduction to the boundary conditions and more detailed code and configurations of this simulation. Step-by-step instruction for building readers’ own simulation is also included.

We thank Yixuan Li and Wern Ng for both carefully reviewing the manuscript prior to submission. This work was supported by the U.K. Engineering and Physical Sciences Research Council through Grant Nos. EP/K037390/1 and EP/M020398/1.

The authors have no conflicts to disclose.

Xiaotian Xu: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Resources (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Michael Newns: Investigation (equal); Software (equal); Validation (supporting); Visualization (supporting); Writing – review & editing (equal). Mark Oxborrow: Funding acquisition (lead); Methodology (equal); Project administration (lead); Resources (lead); Software (lead); Supervision (lead); Writing – review & editing (lead).

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

In COMSOL, each boundary requires a condition (or else several consistent conditions) applied to it. An optimal choice of boundary condition can, furthermore, increase accuracy and reduce the time of calculation.

This section introduces all the choices we have explored so far in our method, which are based on the various “nodes” provided in COMSOL’s weak form PDE interface.

On boundaries that waves cannot (or barely) reach, the fixed boundary condition is a simple choice.

This boundary condition is realized using the Dirichlet Boundary Condition node, where all three input boxes are filled with zeros.

When analyzing some models like an infinitely large substrate, the periodic boundary condition is inevitable. This boundary condition truncates the propagating distance to an integer multiple of the wavelength. It is applied on the boundaries normal to the propagating direction. If one wishes to restrict a simulation to waves that are uniform across a transverse direction, one can implement this restriction by imposing periodic boundary conditions across an analyzed domain that is thin (sub propagation wavelength) in this same direction.

Periodic Condition node is provided by the weak form PDE interface.

On the interface between a piezoelectric material and the environment, the propagation of an acoustic wave stops but the propagation of an electromagnetic wave does not. Thus, the free boundary condition should be applied here for acoustic waves and the continuous boundary condition for the electromagnetic field.

For electromagnetic fields on both sides, the tangential part of $E_$ and the normal part (relative to an interface) of $D_$ are the same,

The free boundary condition can be realized using the default Zero Flux node. The continuity boundary condition can be realized using the Pointwise Constraint node, and its configuration can be found in the supplementary material.

In some models, the piezoelectric device is placed in a cavity enclosed by highly conductive walls. The electric-wall boundary condition can accurately represent such walls while avoiding the need to introduce of a specific material and analyze an additional domain. This condition constrains the tangential part of electric field strength to be zero,

These constraints are applied using three Pointwise Constraint nodes.

An absorbing boundary condition (ABC) is applied to truncate the calculation domain to a reasonable size. Without the ABC, the radiating electromagnetic wave will reflect off the truncating surfaces and disturb (interfere with) wave patterns in other domains. The ABC can reduce this reflection. The ABC is in the form of a constraint and does not need the introduction of an extra cladding domain. The explicit form of the first-order ABC is¹⁶ 

where $r^$ is the spherical radius, $j$ is the imaginary unit, $k$ is the wavenumber, and $s$ is a factor and $s=0.5$ gives the minimum reflection.¹⁶ 

In COMSOL, (A3) is separated into three component directions and written into three Pointwise Constraint nodes.