PDF

Thickness Shear Mode Quartz Oscillator
Introduction
AT cut quartz crystals are widely employed in a range of applications, from oscillators to microbalances. One of the important properties of the AT cut is that the resonant frequency of the crystal is temperature independent to first order. This is desirable in both mass sensing and timing applications. AT cut crystals vibrate in the thickness shear mode — an applied voltage across the faces of the cut produces shear stresses inside the crystal. This example considers the vibration of an AT cut thickness shear oscillator, focusing on the mechanical response of the system in the frequency domain. The effect of a series capacitor on the mechanical resonance is also considered. Adding a series capacitance is a technique frequently employed to tune crystal oscillators.
Model Definition
The model geometry is shown in Figure 1. The oscillator consists of a single (left-handed) quartz disc, supported so as not to impede the motion of the vibrational mode. There are two electrodes on the top and bottom surfaces of the geometry, one of which is grounded.
Figure 1: Model geometry.
In the first version of the model an AC voltage is applied to the top electrode. In the second version, the crystal is still driven by an AC voltage, but a capacitor is placed between the voltage source and the top electrode of the crystal.
Domain Level Equations
Within a piezoelectric crystal there is a coupling between the strain and the electric field, which is determined by the constitutive relation
(1)
Here, S is the strain, T is the stress, E is the electric field, and D is the electric displacement field. The material parameters cE, e, and εS, correspond to the material stiffness, the coupling properties, and the permittivity. These quantities are tensors of rank 4, 3, and 2, respectively, but, since the tensors are highly symmetric for physical reasons, they can be represented as matrices within an abbreviated subscript notation, which is usually more convenient. Equation 1 is implemented by the Piezoelectric Effect branch located under the Multiphysics branch. This constitutive relation is used to couple the equations of Solid Mechanics and Electrostatics, which are solved within the material.
Material Orientation
The orientation of a piezoelectric crystal cut is frequently defined by the system introduced by the IRE standard of 1949 (Ref. 1). This standard has undergone a number of subsequent revisions, with the final revision being the IEEE standard of 1987 (Ref. 2). Unfortunately the 1987 standard contained a number of serious errors and the IEEE subsequently withdrew it. COMSOL Multiphysics therefore adopts the preceding 1978 standard (Ref. 3), which is similar to the 1987 standard, for material property definitions. Most of the material properties in the material library are based on the values given in the book by Auld (Ref. 4), which uses the 1978 IEEE conventions. This is consistent with general practice except in the specific case of quartz, where it is more common to use the 1949 IRE standard to define the material properties. COMSOL Multiphysics therefore provides an additional set of material properties consistent with the 1949 standard for the case of quartz. Note that the material properties for quartz are based on Ref. 5, which uses the 1949 IRE standard (the properties are appropriately modified according to the different standards).
In COMSOL Multiphysics, the stiffness, compliance, coupling, and dielectric material property matrices are defined with respect to the crystal axes. Because the crystal axes for quartz were defined differently between the IRE 1949 and the IEEE 1978 standards (see Figure 3), the signs of several matrix components differ between the two standards (see Table 1). In the absence of a user-defined coordinate system, the crystal axes coincide with the global Xg, Yg, and Zg coordinate axes (the axes shown in the COMSOL Desktop geometry window). When an alternative coordinate system is selected, this system defines the orientation of the crystal axes. This is the mechanism used in COMSOL Multiphysics to define a particular crystal cut, and typically it is necessary to calculate the appropriate Euler angles for the cut (given the thickness orientation for the wafer).
All piezoelectric material properties are defined using the Voigt form of the abbreviated subscript notation, which is universally employed in the literature (this differs from the standard notation used for Linear Elastic Material in the Solid Mechanics interface). The material properties are defined in the material frame, so that if the solid rotates during deformation the material properties rotate with the solid.
Crystal cuts are usually defined by a mechanism introduced by the IEEE/IRE standards. Both standards use a notation that defines the orientation of a virtual slice (the plate) through the crystal. The crystal axes are denoted X, Y, and Z and the plate, which is usually rectangular, is defined as having sides l, w, and t (length, width, and thickness). Initially the plate is aligned with respect to the crystal axes and then up to three rotations are defined, using a right-handed convention about axes embedded along the l, w, and t sides of the plate.
This tutorial uses AT cut quartz, defined in the IEEE 1978 standard as: (YXl-35.25°. The first two letters in the bracketed expression always refer to the initial orientation of the thickness (t) and the length (l) of the plate. Subsequent bracketed letters then define up to three rotational axes, which move with the plate as it is rotated. Angles of rotation about these axes are specified after the bracketed expression in the order of the letters, using a right-handed convention. For AT cut quartz only one rotation, about the l axis, is required. This is illustrated in Figure 2.
Figure 3 shows the final orientation of the plate with respect to the crystal for the AT cut quartz, according to both standards. In this figure, for clarity, the crystal axes labels X, Y, and Z are appended with a subscript “cr” (Xcr, Ycr, and Zcr). Note that within the 1949 IRE Standard, AT cut quartz is defined as: (YXl) +35.25°. As shown in Figure 3, the final orientation of the plate (l, w, and t axes) are not the same between the two standards. Table 2 summarizes the differences between the standards for the AT cut.
When defining the material properties of quartz, the orientation of the crystal axes X, Y, and Z with respect to the crystal differs between the 1978 IEEE standard and the 1949 IRE standard, as shown in Figure 3. (In the figure the crystal axes X, Y, and Z are labeled as Xcr, Ycr, and Zcr for clarity.) A consequence of this is that both the material property matrices and the crystal cut differ between the two standards. Table 1 summarizes the signs for the important matrix elements under the two conventions. Table 2 shows the different definitions of the crystal cuts under the two conventions.
s14
c14
d11
d14
e11
e14
(YXl) +35.25°
(ZXZ: 0°,−35.25°,0°)
(ZXZ: 0°, +54.75°, 0°)
(YXl) -35.25°
(ZXZ: 0°, +35.25°,0°)
(ZXZ: 0°, +125.25°,0°)
Figure 2: Definition of the AT cut of quartz within the IEEE 1978 standard. It is defined as: (YXl) -35.25°. The first two bracketed letters specify the initial orientation of the plate, with the thickness direction (t) along the crystal Y axis, and the length direction (l) along the X axis. Then up to three rotations about axes that move with the plate are specified by the corresponding bracketed letters and the subsequent angles. In this case only one rotation is required about the l axis, of -35.25° (in a right-handed sense).
Figure 3: Summary of two standards for AT cut quartz and the coordinate systems involved in the model setup. The cut plate is in its final orientation (l, w, t), after a single rotation of 35.25° from the crystal axes (Xcr , Ycr , Zcr), according to the standards. The global coordinate system (Xg , Yg , Zg) is given by the orientation of the plate in the model geometry, which can be quite arbitrary as illustrated in the figure.
In general, when setting up a COMSOL model involving crystal cuts, there are three coordinate systems and three sets of relationships among them to be considered, as summarized in Figure 3.
1
2
When creating a COMSOL Multiphysics model, we draw the plate in the geometry window, thereby implicitly defines the relationship between the plate axes (l, w, t) and the global axes (Xg, Yg, Zg). If the plate is oriented differently in the model geometry, then this relationship will also be different.
3
Using this final relationship, we create a local coordinate system that coincides with the crystal axes (Xcr, Ycr, and Zcr), to specify the material orientation, so that the crystal cut in the model is consistent with the standard shown in the figure.
For a specific example, consider the IEEE 1978 standard for AT cut quartz as shown in Figure 2 and Figure 3. As discussed in 2 above, the definition of the appropriate local coordinate system depends on the desired final orientation of the plate in the global coordinate system (how the plate is drawn in the geometry window). One way to set up the plate is to orient its normal (t) along the Yg axis in the global coordinate system, as shown in Figure 4(b). Another way is to orient the plate normal (t) along the global Zg axis, which is the case for the crystal in this model (as in Figure 4(a)). By comparing Figure 4(a) and (b), it is clear that the orientation of the local system (the crystal axes Xcr, Ycr, and Zcr) with respect to the global system (Xg, Yg, and Zg) is different between the two configurations. Therefore it is critical to keep track of this relationship between the local and global system when setting up a model.
There are several methods available to define the local coordinate system with respect to the global system in a model. Usually it is most convenient to define the local coordinates with a Rotated System node, which defines three Euler angles according to the ZXZ convention (rotation about Z, then rotated-X, then rotated-Z). Note that these Euler angles define the rotation starting from the global axes and ending at the local (crystal) axes, as shown in Figure 5. This is distinct from the IRE/IEEE standards that starts the rotation from the crystal (local) axes and ends at the plate axes, as shown in Figure 2.
Figure 4: For the same AT cut of quartz within the IEEE 1978 standard, there are many different possible ways to orient the plate in the model geometry. The thickness direction (t) can be aligned with the global Zg axis (a), or the global Yg axis (b). This results in different Euler angles: 125.25° for the former and 35.25° for the latter. The insets show the plate orientation as seen in the model geometry window.
Figure 5: According to the ZXZ convention adopted by COMSOL Multiphysics, the Euler angles (α, β, γ) define three successive rotations starting from the global axes (Xg, Yg, and Zg) and ending at the local (crystal) axes (Xcr, Ycr, and Zcr).
From Figure 4 it can be seen that for the AT cut quartz, if the plate thickness direction (t) is along the global Zg axis (a), or the global Yg axis (b), then only one non-trivial Euler rotation about the Xg axis is required to bring the global axes (Xg, Yg, and Zg) to align with the crystal (local) axes (Xcr, Ycr, and Zcr). For these two relatively simple cases, the Euler angles can be easily determined: (ZXZ: 0°, +125.25°,0°) for (a) and (ZXZ: 0°, +35.25°,0°) for (b), as shown in the figure and Table 2.
Electrical Circuit
In the first part of the model an AC voltage is applied directly to the top plate of the oscillator, which is grounded. In the second part of the model, a capacitor is added between the voltage source and the oscillator, as shown in Figure 6. In COMSOL Multiphysics, the oscillator is coupled into the circuit using the External I Terminal feature. The terminal boundary condition within the model is set to Circuit and this feature then captures the charge generated by the circuit.
Figure 6: Upper left: Electrical circuit for the first part of the model. Upper right: Electrical circuit for the second part of the model. Bottom: Circuit for the second part of the model as implemented in COMSOL Multiphysics.
Results and Discussion
Figure 7 shows the crystal displacement at its resonant frequency of 5.11 MHz. The form of the displacement shows clearly the shear nature of the resonance. The potential on cut slices through the plate is illustrated in Figure 8. The mechanical domain frequency response of the oscillator is shown in Figure 9. A clear resonance is apparent, with a resonant frequency close to 5.11 MHz.
The addition of a series capacitance between the oscillator and the voltage source is expected to pull the resonant frequency to higher values. Figure 10 shows this effect, with the resonant frequency increasing the most for smaller values of the series capacitance (in the limit of very large series capacitance the impedance of the series capacitor goes to zero and produces the result shown in Figure 9).
Figure 7: Displacement of the crystal at resonance.
Figure 8: Electric potential inside the crystal at resonance.
Figure 9: Mechanical response of the structure with no series capacitance.
Figure 10: Mechanical response of the structure with different series capacitances.
References
1. “Standards on Piezoelectric Crystals, 1949”, Proceedings of the I. R. E.,vol. 37, no. 12, pp. 1378–1395, 1949.
2. IEEE Standard on Piezoelectricity, ANSI/IEEE Standard 176-1987, 1987.
3. IEEE Standard on Piezoelectricity, ANSI/IEEE Standard 176-1978, 1978.
4. B. A. Auld, Acoustic Fields and Waves in Solids, Krieger Publishing Company, 1990.
5. R. Bechmann, “Elastic and Piezoelectric Constants of Alpha-Quartz”, Physical Review B, vol. 110 no. 5, pp. 1060–1061, 1958.
Application Library path: MEMS_Module/Piezoelectric_Devices/thickness_shear_quartz_oscillator
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  3D.
2
In the Select Physics tree, select Structural Mechanics>Electromagnetics-Structure Interaction>Piezoelectricity>Piezoelectricity, Solid.
3
Click Add.
We will select the Adaptive Frequency Sweep study step later.
4
Geometry 1
Add parameters for the model geometry and series capacitance.
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Instead of building the geometry from scratch, import a pre-built mesh. This showcases a less-known functionality of COMSOL Multiphysics to create and solve a model from imported mesh without needing any geometry object to be created.
To keep the computational time and RAM requirements as low as possible, we use a coarse mesh. For this kind of high-Q systems, a much finer mesh is required for the resonant structure to be truly converged.
Mesh 1
Import 1
1
In the Mesh toolbar, click  Import.
2
In the Settings window for Import, locate the Import section.
3
Click Browse.
4
5
Click Import.
Create an operator to evaluate quantities on a point for the Adaptive Frequency Sweep study step.
Definitions
Integration 1 (intop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Integration.
2
In the Settings window for Integration, locate the Source Selection section.
3
From the Geometric entity level list, choose Point.
4
The Adaptive Frequency Sweep study step will generate a high resolution frequency sweep. To avoid large file size, create an "explicit selection" to store solution data only on the external surfaces of the modeling domain.
All surfaces
1
In the Definitions toolbar, click  Explicit.
2
In the Settings window for Explicit, type All surfaces in the Label text field.
3
Locate the Input Entities section. From the Geometric entity level list, choose Boundary.
4
Select the All boundaries check box.
Set up a rotated system appropriate for AT cut Quartz.
Rotated System 2 (sys2)
1
In the Definitions toolbar, click  Coordinate Systems and choose Rotated System.
2
In the Settings window for Rotated System, locate the Rotation section.
3
Find the Euler angles (Z-X-Z) subsection. In the β text field, type 125.25[deg].
Add Material
1
In the Home toolbar, click  Add Material to open the Add Material window.
2
Go to the Add Material window.
3
In the tree, select Piezoelectric>Quartz LH (1978 IEEE).
4
Click Add to Component in the window toolbar.
5
In the Home toolbar, click  Add Material to close the Add Material window.
Solid Mechanics (solid)
Use the rotated system to define the orientation of the crystal.
Piezoelectric Material 1
1
In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid) click Piezoelectric Material 1.
2
In the Settings window for Piezoelectric Material, locate the Coordinate System Selection section.
3
From the Coordinate system list, choose Rotated System 2 (sys2).
Add damping to the model.
Mechanical Damping 1
1
In the Physics toolbar, click  Attributes and choose Mechanical Damping.
2
In the Settings window for Mechanical Damping, locate the Damping Settings section.
3
From the Damping type list, choose Isotropic loss factor.
4
From the ηs list, choose User defined. In the associated text field, type 1e-3.
Electrostatics (es)
Add electrical boundary conditions to the model. First add a Terminal boundary condition that connects the electrode to an external circuit.
1
In the Model Builder window, under Component 1 (comp1) click Electrostatics (es).
Terminal 1
1
In the Physics toolbar, click  Boundaries and choose Terminal.
2
3
In the Settings window for Terminal, locate the Terminal section.
4
From the Terminal type list, choose Circuit.
Override the preceding boundary condition with a constant potential boundary condition to compute the response of the device without a series capacitance. This node will be disabled in the study when the circuit is included in the model.
Terminal 2
1
In the Physics toolbar, click  Boundaries and choose Terminal.
2
In the Settings window for Terminal, locate the Terminal section.
3
In the Terminal name text field, type 1.
4
5
From the Terminal type list, choose Voltage.
6
In the V0 text field, type 10.
Ground 1
1
In the Physics toolbar, click  Boundaries and choose Ground.
2
Root
In the second study in this model, the effect of a series capacitor on the device response will be investigated. Add an Electrical Circuit interface to model the capacitor.
Add Physics
1
In the Physics toolbar, click  Add Physics to open the Add Physics window.
2
Go to the Add Physics window.
3
In the tree, select AC/DC>Electrical Circuit (cir).
4
Click Add to Component 1 in the window toolbar.
5
In the Home toolbar, click  Add Physics to close the Add Physics window.
Electrical Circuit (cir)
Features in the Electric Circuits interface are connected by specifying connecting node numbers for each port of the device.
Ground Node 1 (gnd1)
A ground node is automatically added to the circuit, with the default node number of 0.
Next add a voltage source between the ground node and a (newly created) node with number 2.
1
In the Model Builder window, click Electrical Circuit (cir).
Voltage Source 1 (V1)
1
In the Electrical Circuit toolbar, click  Voltage Source.
2
In the Settings window for Voltage Source, locate the Node Connections section.
3
4
Locate the Device Parameters section. From the Source type list, choose AC-source.
5
In the vsrc text field, type 10.
Add a capacitor between the voltage source output (node 2) and a new node, 1.
Capacitor 1 (C1)
1
In the Electrical Circuit toolbar, click  Capacitor.
2
In the Settings window for Capacitor, locate the Node Connections section.
3
4
Locate the Device Parameters section. In the C text field, type Cs.
Connect node 1 to the Terminal feature in the model using the External I-Terminal feature.
External I-terminal 1 (termI1)
1
In the Electrical Circuit toolbar, click  External I-terminal.
Couple the electric potential from the Terminal in the electrostatics interface back into the model.
2
In the Settings window for External I-terminal, locate the Node Connections section.
3
In the Node name text field, type 1.
4
Locate the External Terminal section. From the V list, choose Terminal voltage (es).
Mesh 1
The mesh used here is somewhat coarse for this kind of high-Q systems. The reason is mainly to keep the computational time and RAM requirements as low as possible. Interested users are encouraged to solve the problem for finer mesh settings.
Set up and solve an Adaptive Frequency Sweep study, which is optimized for resolving narrow resonant peaks without excessive computation.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select Empty Study.
4
Click Add Study in the window toolbar.
5
In the Home toolbar, click  Add Study to close the Add Study window.
Study 1
Adaptive Frequency Sweep
1
In the Study toolbar, click  Study Steps and choose Frequency Domain>Adaptive Frequency Sweep.
2
In the Settings window for Adaptive Frequency Sweep, locate the Study Settings section.
3
In the Frequencies text field, type range(5.095[MHz],0.2[kHz],5.13[MHz]).
4
From the AWE expression type list, choose User controlled.
5
In the first study, disable the electrical circuit.
6
Locate the Physics and Variables Selection section. In the table, clear the Solve for check box for Electrical Circuit (cir).
7
In the Study toolbar, click  Compute.
Results
Instead of the stress plot, visualize the mode shape of the device at resonance.
Displacement
1
In the Model Builder window, under Results click Stress (solid).
2
In the Settings window for 3D Plot Group, type Displacement in the Label text field.
3
Locate the Data section. From the Parameter value (freq (Hz)) list, choose 5.112E6.
Surface 1
1
In the Model Builder window, expand the Displacement node, then click Surface 1.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Solid Mechanics>Displacement>solid.disp - Displacement magnitude - m.
3
Locate the Expression section. From the Unit list, choose nm.
4
Click the  Zoom Extents button in the Graphics toolbar.
5
In the Displacement toolbar, click  Plot.
The second default plot shows the electric potential within the device. For a better view, plot 5 slices in xy planes.
Electric Potential (es)
1
In the Model Builder window, click Electric Potential (es).
2
In the Settings window for 3D Plot Group, locate the Data section.
3
From the Parameter value (freq (Hz)) list, choose 5.112E6.
Multislice 1
1
In the Model Builder window, expand the Electric Potential (es) node, then click Multislice 1.
2
In the Settings window for Multislice, locate the Multiplane Data section.
3
Find the x-planes subsection. In the Planes text field, type 5.
4
Find the y-planes subsection. In the Planes text field, type 0.
5
Find the z-planes subsection. In the Planes text field, type 0.
6
In the Electric Potential (es) toolbar, click  Plot.
Add a plot to show the mechanical response of the device.
Mechanical Response
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, type Mechanical Response in the Label text field.
Point Graph 1
1
Right-click Mechanical Response and choose Point Graph.
2
3
In the Settings window for Point Graph, locate the y-Axis Data section.
4
In the Expression text field, type abs(u).
5
From the Unit list, choose nm.
6
Locate the x-Axis Data section. From the Unit list, choose MHz.
7
In the Mechanical Response toolbar, click  Plot.
Now set up a study to compute the frequency response of the device with different capacitors added in series.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select Empty Study.
4
Click Add Study in the window toolbar.
5
In the Home toolbar, click  Add Study to close the Add Study window.
Study 2
Adaptive Frequency Sweep
1
In the Study toolbar, click  Study Steps and choose Frequency Domain>Adaptive Frequency Sweep.
2
In the Settings window for Adaptive Frequency Sweep, locate the Study Settings section.
3
In the Frequencies text field, type range(5.095[MHz],0.2[kHz],5.13[MHz]).
4
From the AWE expression type list, choose User controlled.
5
6
Locate the Physics and Variables Selection section. Select the Modify model configuration for study step check box.
7
In the Physics and variables selection tree, select Component 1 (comp1)>Electrostatics (es)>Terminal 2.
8
Click  Disable.
To reduce file size, only store solution data on the external surfaces.
9
Locate the Values of Dependent Variables section. Find the Store fields in output subsection. From the Settings list, choose For selections.
10
Under Selections, click  Add.
11
In the Add dialog box, select All surfaces in the Selections list.
12
Parametric Sweep
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
4
5
In the Model Builder window, click Study 2.
6
In the Settings window for Study, locate the Study Settings section.
7
Clear the Generate default plots check box.
8
In the Study toolbar, click  Compute.
Results
Re-plot the mechanical response with the additional series capacitance.
Mechanical response, Parametric
1
In the Model Builder window, right-click Mechanical Response and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Mechanical response, Parametric in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 2/Parametric Solutions 1 (sol3).
4
Locate the Legend section. From the Position list, choose Upper left.
Point Graph 1
1
In the Model Builder window, expand the Mechanical response, Parametric node, then click Point Graph 1.
2
In the Settings window for Point Graph, click to expand the Legends section.
3
Select the Show legends check box.
4
In the Mechanical response, Parametric toolbar, click  Plot.
Note how the mechanical resonant frequency is ’pulled’ by the series capacitance.