Appendix 2: Film Bulk Acoustic Resonator (FBAR)

Generated from FBAR.wxm with rfMaxima 0.2.5 on Maxima 5.46.0. All results were computed by executing the notebook; none is transcribed.

(See: Microsystem Design by Stephen D. Senturia, pp. 192-193, 570-577)
(See: An Introduction to the Theory of Piezoelectricity by Jiashi Yang, pp. 89-93)

Linear Piezoelectricity

Displacement (u), Strain (S), and Stress (T)

u: displacement [m]

u:matrix([u_x(x,y,z,t)],[u_y(x,y,z,t)],[u_z(x,y,z,t)])$
u2S(u):=matrix([diff(u[1,1],x)],[diff(u[2,1],y)],[diff(u[3,1],z)],[diff(u[2,1],z)+diff(u[3,1],y)],[diff(u[1,1],z)+diff(u[3,1],x)],[diff(u[1,1],y)+diff(u[2,1],x)])$

S: strain [] (%epsilon: normal strain, %gamma: shear strain)

S:matrix([%epsilon_x(x,y,z,t)],[%epsilon_y(x,y,z,t)],[%epsilon_z(x,y,z,t)],[%gamma_yz(x,y,z,t)],[%gamma_zx(x,y,z,t)],[%gamma_xy(x,y,z,t)])$
S=u2S(u);

\begin{pmatrix} \epsilon_{x}\left(x , y , z , t\right)\\ \epsilon_{y}\left(x , y , z , t\right)\\ \epsilon_{z}\left(x , y , z , t\right)\\ \gamma_{yz}\left(x , y , z , t\right)\\ \gamma_{zx}\left(x , y , z , t\right)\\ \gamma_{xy}\left(x , y , z , t\right) \end{pmatrix} =\begin{pmatrix} {{d}\over{d\,x}}\,\mathit{u}_{x}\left(x , y , z , t\right)\\ {{d}\over{d\,y}}\,\mathit{u}_{y}\left(x , y , z , t\right)\\ {{d}\over{d\,z}}\,\mathit{u}_{z}\left(x , y , z , t\right)\\ {{d}\over{d\,y}}\,\mathit{u}_{z}\left(x , y , z , t\right)+{{d}\over{d\,z}}\,\mathit{u}_{y}\left(x , y , z , t\right)\\ {{d}\over{d\,x}}\,\mathit{u}_{z}\left(x , y , z , t\right)+{{d}\over{d\,z}}\,\mathit{u}_{x}\left(x , y , z , t\right)\\ {{d}\over{d\,x}}\,\mathit{u}_{y}\left(x , y , z , t\right)+{{d}\over{d\,y}}\,\mathit{u}_{x}\left(x , y , z , t\right) \end{pmatrix}

grad(u):=matrix([diff(u[1,1],x,1),diff(u[1,1],y,1),diff(u[1,1],z,1)],[diff(u[2,1],x,1),diff(u[2,1],y,1),diff(u[2,1],z,1)],[diff(u[3,1],x,1),diff(u[3,1],y,1),diff(u[3,1],z,1)])$
grad(u);
1/2*(grad(u)+transpose(grad(u)));

\begin{pmatrix} {{d}\over{d\,x}}\,\mathit{u}_{x}\left(x , y , z , t\right)&{{d}\over{d\,y}}\,\mathit{u}_{x}\left(x , y , z , t\right)&{{d}\over{d\,z}}\,\mathit{u}_{x}\left(x , y , z , t\right)\\ {{d}\over{d\,x}}\,\mathit{u}_{y}\left(x , y , z , t\right)&{{d}\over{d\,y}}\,\mathit{u}_{y}\left(x , y , z , t\right)&{{d}\over{d\,z}}\,\mathit{u}_{y}\left(x , y , z , t\right)\\ {{d}\over{d\,x}}\,\mathit{u}_{z}\left(x , y , z , t\right)&{{d}\over{d\,y}}\,\mathit{u}_{z}\left(x , y , z , t\right)&{{d}\over{d\,z}}\,\mathit{u}_{z}\left(x , y , z , t\right) \end{pmatrix}

\begin{pmatrix} {{d}\over{d\,x}}\,\mathit{u}_{x}\left(x , y , z , t\right)&{{{{d}\over{d\,x}}\,\mathit{u}_{y}\left(x , y , z , t\right)+{{d}\over{d\,y}}\,\mathit{u}_{x}\left(x , y , z , t\right)}\over{2}}&{{{{d}\over{d\,x}}\,\mathit{u}_{z}\left(x , y , z , t\right)+{{d}\over{d\,z}}\,\mathit{u}_{x}\left(x , y , z , t\right)}\over{2}}\\ {{{{d}\over{d\,x}}\,\mathit{u}_{y}\left(x , y , z , t\right)+{{d}\over{d\,y}}\,\mathit{u}_{x}\left(x , y , z , t\right)}\over{2}}&{{d}\over{d\,y}}\,\mathit{u}_{y}\left(x , y , z , t\right)&{{{{d}\over{d\,y}}\,\mathit{u}_{z}\left(x , y , z , t\right)+{{d}\over{d\,z}}\,\mathit{u}_{y}\left(x , y , z , t\right)}\over{2}}\\ {{{{d}\over{d\,x}}\,\mathit{u}_{z}\left(x , y , z , t\right)+{{d}\over{d\,z}}\,\mathit{u}_{x}\left(x , y , z , t\right)}\over{2}}&{{{{d}\over{d\,y}}\,\mathit{u}_{z}\left(x , y , z , t\right)+{{d}\over{d\,z}}\,\mathit{u}_{y}\left(x , y , z , t\right)}\over{2}}&{{d}\over{d\,z}}\,\mathit{u}_{z}\left(x , y , z , t\right) \end{pmatrix}

T: stress, force per unit area [N/m^2] (%sigma: normall stress, %tau: shear stress)

T:matrix([%sigma_x(x,y,z,t)],[%sigma_y(x,y,z,t)],[%sigma_z(x,y,z,t)],[%tau_yz(x,y,z,t)],[%tau_zx(x,y,z,t)],[%tau_xy(x,y,z,t)]);
divT(T):=matrix([diff(T[1,1],x,1)+diff(T[6,1],y,1)+diff(T[5,1],z,1)],[diff(T[6,1],x,1)+diff(T[2,1],y,1)+diff(T[4,1],z,1)],[diff(T[5,1],x,1)+diff(T[4,1],y,1)+diff(T[3,1],z,1)])$
divT(T);

\begin{pmatrix} \sigma_{x}\left(x , y , z , t\right)\\ \sigma_{y}\left(x , y , z , t\right)\\ \sigma_{z}\left(x , y , z , t\right)\\ \tau_{yz}\left(x , y , z , t\right)\\ \tau_{zx}\left(x , y , z , t\right)\\ \tau_{xy}\left(x , y , z , t\right) \end{pmatrix}

\begin{pmatrix} {{d}\over{d\,z}}\,\tau_{zx}\left(x , y , z , t\right)+{{d}\over{d\,y}}\,\tau_{xy}\left(x , y , z , t\right)+{{d}\over{d\,x}}\,\sigma_{x}\left(x , y , z , t\right)\\ {{d}\over{d\,z}}\,\tau_{yz}\left(x , y , z , t\right)+{{d}\over{d\,x}}\,\tau_{xy}\left(x , y , z , t\right)+{{d}\over{d\,y}}\,\sigma_{y}\left(x , y , z , t\right)\\ {{d}\over{d\,x}}\,\tau_{zx}\left(x , y , z , t\right)+{{d}\over{d\,y}}\,\tau_{yz}\left(x , y , z , t\right)+{{d}\over{d\,z}}\,\sigma_{z}\left(x , y , z , t\right) \end{pmatrix}

Electric Displacement (D) and Electric Field (E)

D: electric displacement [C/m^2]

D:matrix([D_x(x,y,z,t)],[D_y(x,y,z,t)],[D_z(x,y,z,t)]);
/* Divergence of D is a SCALAR.  0.2.4 returned a 3x1 column of the three
   separate terms, which then made Gauss's law read as three equations. */
divD(D):=diff(D[1,1],x,1)+diff(D[2,1],y,1)+diff(D[3,1],z,1)$
divD(D);

\begin{pmatrix} \mathit{D}_{x}\left(x , y , z , t\right)\\ \mathit{D}_{y}\left(x , y , z , t\right)\\ \mathit{D}_{z}\left(x , y , z , t\right) \end{pmatrix}

{{d}\over{d\,z}}\,\mathit{D}_{z}\left(x , y , z , t\right)+{{d}\over{d\,y}}\,\mathit{D}_{y}\left(x , y , z , t\right)+{{d}\over{d\,x}}\,\mathit{D}_{x}\left(x , y , z , t\right)

E: electric field [V/m]

E:matrix([E_x(x,y,z,t)],[E_y(x,y,z,t)],[E_z(x,y,z,t)]);

\begin{pmatrix} \mathit{E}_{x}\left(x , y , z , t\right)\\ \mathit{E}_{y}\left(x , y , z , t\right)\\ \mathit{E}_{z}\left(x , y , z , t\right) \end{pmatrix}

Linear Constitutive Tensors

The number of independent non-zero constants depends on the crystal symmetry. AlN and ZnO are materials of the 6mm crystal class.

%epsilon: dielectric permittivity tensor [C/(V.m)] (at constant strain)

declare([%epsilon_0,%epsilon_11,%epsilon_33],[constant,real,scalar])$
%epsilon:matrix([%epsilon_11,0,0],[0,%epsilon_11,0],[0,0,%epsilon_33]);

\begin{pmatrix} \epsilon_{11}&0&0\\ 0&\epsilon_{11}&0\\ 0&0&\epsilon_{33} \end{pmatrix}

C: elastic stiffness tensor [N/m^2] (at constant electric field)

declare([C_11,C_12,C_13,C_33,C_44],[constant,real,scalar])$
C:matrix([C_11, C_12, C_13, 0, 0, 0],[C_12, C_11, C_13, 0, 0, 0],[C_13, C_13, C_33, 0, 0, 0],[0, 0, 0, C_44, 0, 0],[0, 0, 0, 0, C_44, 0],[0, 0, 0, 0, 0, (C_11-C_12)/2]);

\begin{pmatrix} C_{11}&C_{12}&C_{13}&0&0&0\\ C_{12}&C_{11}&C_{13}&0&0&0\\ C_{13}&C_{13}&C_{33}&0&0&0\\ 0&0&0&C_{44}&0&0\\ 0&0&0&0&C_{44}&0\\ 0&0&0&0&0&{{C_{11}-C_{12}}\over{2}} \end{pmatrix}

e: piezoelectricity tensor [C/m^2]

declare([e_15,e_31,e_33],[constant,real,scalar])$
e:matrix([0,0,0,0,e_15,0],[0,0,0,e_15,0,0],[e_31,e_31,e_33,0,0,0]);

\begin{pmatrix} 0&0&0&0&e_{15}&0\\ 0&0&0&e_{15}&0&0\\ e_{31}&e_{31}&e_{33}&0&0&0 \end{pmatrix}

Lippmann's Linear Constitutive Equations in Charge-Stress Form

T:C.S-transpose(e).E;
D:e.S+%epsilon.E;

\begin{pmatrix} -e_{31}\,\mathit{E}_{z}\left(x , y , z , t\right)+C_{13}\,\epsilon_{z}\left(x , y , z , t\right)+C_{12}\,\epsilon_{y}\left(x , y , z , t\right)+C_{11}\,\epsilon_{x}\left(x , y , z , t\right)\\ -e_{31}\,\mathit{E}_{z}\left(x , y , z , t\right)+C_{13}\,\epsilon_{z}\left(x , y , z , t\right)+C_{11}\,\epsilon_{y}\left(x , y , z , t\right)+C_{12}\,\epsilon_{x}\left(x , y , z , t\right)\\ -e_{33}\,\mathit{E}_{z}\left(x , y , z , t\right)+C_{33}\,\epsilon_{z}\left(x , y , z , t\right)+C_{13}\,\epsilon_{y}\left(x , y , z , t\right)+C_{13}\,\epsilon_{x}\left(x , y , z , t\right)\\ C_{44}\,\gamma_{yz}\left(x , y , z , t\right)-e_{15}\,\mathit{E}_{y}\left(x , y , z , t\right)\\ C_{44}\,\gamma_{zx}\left(x , y , z , t\right)-e_{15}\,\mathit{E}_{x}\left(x , y , z , t\right)\\ {{\left(C_{11}-C_{12}\right)\,\gamma_{xy}\left(x , y , z , t\right)}\over{2}} \end{pmatrix}

\begin{pmatrix} \epsilon_{11}\,\mathit{E}_{x}\left(x , y , z , t\right)+e_{15}\,\gamma_{zx}\left(x , y , z , t\right)\\ \epsilon_{11}\,\mathit{E}_{y}\left(x , y , z , t\right)+e_{15}\,\gamma_{yz}\left(x , y , z , t\right)\\ \epsilon_{33}\,\mathit{E}_{z}\left(x , y , z , t\right)+e_{33}\,\epsilon_{z}\left(x , y , z , t\right)+e_{31}\,\epsilon_{y}\left(x , y , z , t\right)+e_{31}\,\epsilon_{x}\left(x , y , z , t\right) \end{pmatrix}

Boundary Value Problem

Assume a thickness-stretch vibration of an AlN plate, excited by an electric field along the Z-axis. Assume the plate can not move along the X- and Y-axis: u_x(z,t)=0, u_y(z,t)=0.

%epsilon_x(x,y,z,t):=0$
%epsilon_y(x,y,z,t):=0$
%epsilon_z(x,y,z,t):=diff(u_z(z,t),z,1)$
%gamma_yz(x,y,z,t):=0$
%gamma_zx(x,y,z,t):=0$
%gamma_xy(x,y,z,t):=0$

Assume the X- and Y-component of the electric field to be zero: E_x(x,y,z,t)=0, E_y(x,y,z,t)=0.

E:matrix([0],[0],[E_z(z,t)]);
T:ev(C.S-transpose(e).E);
D:ev(e.S+%epsilon.E);

\begin{pmatrix} 0\\ 0\\ \mathit{E}_{z}\left(z , t\right) \end{pmatrix}

\begin{pmatrix} C_{13}\,\left({{d}\over{d\,z}}\,\mathit{u}_{z}\left(z , t\right)\right)-e_{31}\,\mathit{E}_{z}\left(z , t\right)\\ C_{13}\,\left({{d}\over{d\,z}}\,\mathit{u}_{z}\left(z , t\right)\right)-e_{31}\,\mathit{E}_{z}\left(z , t\right)\\ C_{33}\,\left({{d}\over{d\,z}}\,\mathit{u}_{z}\left(z , t\right)\right)-e_{33}\,\mathit{E}_{z}\left(z , t\right)\\ 0\\ 0\\ 0 \end{pmatrix}

\begin{pmatrix} 0\\ 0\\ e_{33}\,\left({{d}\over{d\,z}}\,\mathit{u}_{z}\left(z , t\right)\right)+\epsilon_{33}\,\mathit{E}_{z}\left(z , t\right) \end{pmatrix}

Potential Energy

Electrostatic Energy Density, w_e

w_e:expand(ev(1/2*transpose(E).D));

{{e_{33}\,\mathit{E}_{z}\left(z , t\right)\,\left({{d}\over{d\,z}}\,\mathit{u}_{z}\left(z , t\right)\right)}\over{2}}+{{\epsilon_{33}\,\mathit{E}_{z}\left(z , t\right)^2}\over{2}}

Elastic Energy Density, w_m

w_m:expand(ev(1/2*transpose(S).C.S));

{{C_{33}\,\left({{d}\over{d\,z}}\,\mathit{u}_{z}\left(z , t\right)\right)^2}\over{2}}

Total Energy Density

w_t:w_e+w_m;

{{C_{33}\,\left({{d}\over{d\,z}}\,\mathit{u}_{z}\left(z , t\right)\right)^2}\over{2}}+{{e_{33}\,\mathit{E}_{z}\left(z , t\right)\,\left({{d}\over{d\,z}}\,\mathit{u}_{z}\left(z , t\right)\right)}\over{2}}+{{\epsilon_{33}\,\mathit{E}_{z}\left(z , t\right)^2}\over{2}}

Newton's Law of Motion for a Continuum [PDE in AlN Volume]

Force density

ev(T);
f:ratsimp(divT(ev(T)));

\begin{pmatrix} C_{13}\,\left({{d}\over{d\,z}}\,\mathit{u}_{z}\left(z , t\right)\right)-e_{31}\,\mathit{E}_{z}\left(z , t\right)\\ C_{13}\,\left({{d}\over{d\,z}}\,\mathit{u}_{z}\left(z , t\right)\right)-e_{31}\,\mathit{E}_{z}\left(z , t\right)\\ C_{33}\,\left({{d}\over{d\,z}}\,\mathit{u}_{z}\left(z , t\right)\right)-e_{33}\,\mathit{E}_{z}\left(z , t\right)\\ 0\\ 0\\ 0 \end{pmatrix}

\begin{pmatrix} 0\\ 0\\ C_{33}\,\left({{d^2}\over{d\,z^2}}\,\mathit{u}_{z}\left(z , t\right)\right)-e_{33}\,\left({{d}\over{d\,z}}\,\mathit{E}_{z}\left(z , t\right)\right) \end{pmatrix}

Viscosity and Q Factor

%eta: viscosity tensor [N/(m^2.s)] %tau_m is the motional time constant %omega_m is the mechanical resonant radian frequency Q_m is the mechanical Q factor

%tau_m:1/(%omega_m*Q_m);
%eta:%tau_m*C;

{{1}\over{\omega_{m}\,\mathit{Q}_{m}}}

\begin{pmatrix} {{C_{11}}\over{\omega_{m}\,\mathit{Q}_{m}}}&{{C_{12}}\over{\omega_{m}\,\mathit{Q}_{m}}}&{{C_{13}}\over{\omega_{m}\,\mathit{Q}_{m}}}&0&0&0\\ {{C_{12}}\over{\omega_{m}\,\mathit{Q}_{m}}}&{{C_{11}}\over{\omega_{m}\,\mathit{Q}_{m}}}&{{C_{13}}\over{\omega_{m}\,\mathit{Q}_{m}}}&0&0&0\\ {{C_{13}}\over{\omega_{m}\,\mathit{Q}_{m}}}&{{C_{13}}\over{\omega_{m}\,\mathit{Q}_{m}}}&{{C_{33}}\over{\omega_{m}\,\mathit{Q}_{m}}}&0&0&0\\ 0&0&0&{{C_{44}}\over{\omega_{m}\,\mathit{Q}_{m}}}&0&0\\ 0&0&0&0&{{C_{44}}\over{\omega_{m}\,\mathit{Q}_{m}}}&0\\ 0&0&0&0&0&{{C_{11}-C_{12}}\over{2\,\omega_{m}\,\mathit{Q}_{m}}} \end{pmatrix}

Newton's Law of Motion for a Continuum

%rho: mass density [kg/m^3]

f=%rho*diff(u,t,2);

\begin{pmatrix} 0\\ 0\\ C_{33}\,\left({{d^2}\over{d\,z^2}}\,\mathit{u}_{z}\left(z , t\right)\right)-e_{33}\,\left({{d}\over{d\,z}}\,\mathit{E}_{z}\left(z , t\right)\right) \end{pmatrix} =\begin{pmatrix} \rho\,\left({{d^2}\over{d\,t^2}}\,\mathit{u}_{x}\left(x , y , z , t\right)\right)\\ \rho\,\left({{d^2}\over{d\,t^2}}\,\mathit{u}_{y}\left(x , y , z , t\right)\right)\\ \rho\,\left({{d^2}\over{d\,t^2}}\,\mathit{u}_{z}\left(x , y , z , t\right)\right) \end{pmatrix}

f_z:ev(f[3,1]);
f_z=%rho*diff(u_z(z,t),t,2);

C_{33}\,\left({{d^2}\over{d\,z^2}}\,\mathit{u}_{z}\left(z , t\right)\right)-e_{33}\,\left({{d}\over{d\,z}}\,\mathit{E}_{z}\left(z , t\right)\right)

C_{33}\,\left({{d^2}\over{d\,z^2}}\,\mathit{u}_{z}\left(z , t\right)\right)-e_{33}\,\left({{d}\over{d\,z}}\,\mathit{E}_{z}\left(z , t\right)\right)=\rho\,\left({{d^2}\over{d\,t^2}}\,\mathit{u}_{z}\left(z , t\right)\right)

Gauss' Law [PDE in AlN Volume]

divD(ev(D))=0;

e_{33}\,\left({{d^2}\over{d\,z^2}}\,\mathit{u}_{z}\left(z , t\right)\right)+\epsilon_{33}\,\left({{d}\over{d\,z}}\,\mathit{E}_{z}\left(z , t\right)\right)=0

Calculating 'diff(E_z(z,t),z,1)

-e_33*(diff(u_z(z,t),z,2))/%epsilon_33;

-{{e_{33}\,\left({{d^2}\over{d\,z^2}}\,\mathit{u}_{z}\left(z , t\right)\right)}\over{\epsilon_{33}}}

Calculating E_z(z,t) B_1 is an integration constant.

declare([B_1],[constant,real,scalar])$
define(E_z(z,t),-(e_33*(diff(u_z(z,t),z,1)))/%epsilon_33-B_1);

\mathit{E}_{z}\left(z , t\right):=-{{e_{33}\,\left({{d}\over{d\,z}}\,\mathit{u}_{z}\left(z , t\right)\right)}\over{\epsilon_{33}}}-B_{1}

Boundary Conditions at the Plate Surface

Assume a solution u_z(z,t):=A_1*sin(%xi*z)*exp(%i*%omega_m*t) with %xi:sqrt(%rho/(C_33*(1+k_33^2)))*%omega
k_33: sqrt(e_33^2/(%epsilon_33*C_33))

declare([A_1],[constant,real,scalar])$
u_z(z):=A_1*sin(%xi*z);
ratsimp(diff(u_z(z),z,1));
ratsimp(diff(u_z(z),z,2));

\mathit{u}_{z}\left(z\right):=A_{1}\,\sin \left(\xi\,z\right)

A_{1}\,\xi\,\cos \left(\xi\,z\right)

-A_{1}\,\xi^2\,\sin \left(\xi\,z\right)

define(T_z(z),C_33*(1+k_33^2)*%th(2)+e_33*B_1);
D_z(z):=-%epsilon_33*B_1;
define(E_z(z),-(e_33*%th(4))/%epsilon_33-B_1);
define(f_z(z),C_33*(1+k_33^2)*%th(4));
define(f_s_z(z),C_33*%th(5));
define(f_p_z(z),-e_33*(diff(E_z(z),z,1)));

\mathit{T}_{z}\left(z\right):=A_{1}\,C_{33}\,\xi\,\left(k_{33}^2+1\right)\,\cos \left(\xi\,z\right)+B_{1}\,e_{33}

\mathit{D}_{z}\left(z\right):=\left(-\epsilon_{33}\right)\,B_{1}

\mathit{E}_{z}\left(z\right):=-{{A_{1}\,e_{33}\,\xi\,\cos \left(\xi\,z\right)}\over{\epsilon_{33}}}-B_{1}

\mathit{f}_{z}\left(z\right):=-A_{1}\,C_{33}\,\xi^2\,\left(k_{33}^2+1\right)\,\sin \left(\xi\,z\right)

\mathit{f}_{s,z}\left(z\right):=-A_{1}\,C_{33}\,\xi^2\,\sin \left(\xi\,z\right)

\mathit{f}_{p,z}\left(z\right):=-{{A_{1}\,e_{33}^2\,\xi^2\,\sin \left(\xi\,z\right)}\over{\epsilon_{33}}}

Calculating the potential %Phi(z) B_2 is an integration constant.

declare([B_2],[constant,real,scalar])$
%Phi(z):=e_33*u_z(z)/%epsilon_33+B_1*z+B_2;

\Phi\left(z\right):={{e_{33}\,\mathit{u}_{z}\left(z\right)}\over{\epsilon_{33}}}+B_{1}\,z+B_{2}

%sigma_z(z=-g_0/2)=0, %sigma_z(z=g_0/2)=0, %Phi(z=g_0/2)-%Phi(z=-g_0/2)=V with, g_0: thickness [m] V: applied voltage [V]

solve([ev(T_z(g_0/2))=0,ev(%Phi(g_0/2)-%Phi(-g_0/2))=V],[A_1,B_1])$
A_1:ratsimp(ev(A_1,%));
B_1:ratsimp(ev(B_1,%th(2)));
ratsimp(ev(A_1,%xi=%pi/g_0));
ratsimp(ev(B_1,%xi=%pi/g_0));

-{{\epsilon_{33}\,e_{33}\,V}\over{\epsilon_{33}\,C_{33}\,\xi\,g_{0}\,\cos \left({{\xi\,g_{0}}\over{2}}\right)\,k_{33}^2-2\,e_{33}^2\,\sin \left({{\xi\,g_{0}}\over{2}}\right)+\epsilon_{33}\,C_{33}\,\xi\,g_{0}\,\cos \left({{\xi\,g_{0}}\over{2}}\right)}}

{{\epsilon_{33}\,C_{33}\,\xi\,V\,\cos \left({{\xi\,g_{0}}\over{2}}\right)\,k_{33}^2+\epsilon_{33}\,C_{33}\,\xi\,V\,\cos \left({{\xi\,g_{0}}\over{2}}\right)}\over{\epsilon_{33}\,C_{33}\,\xi\,g_{0}\,\cos \left({{\xi\,g_{0}}\over{2}}\right)\,k_{33}^2-2\,e_{33}^2\,\sin \left({{\xi\,g_{0}}\over{2}}\right)+\epsilon_{33}\,C_{33}\,\xi\,g_{0}\,\cos \left({{\xi\,g_{0}}\over{2}}\right)}}

{{\epsilon_{33}\,V}\over{2\,e_{33}}}

0

Elastic Wave Equation and Christoffel Equation

(See: https://sepwww.stanford.edu/data/media/public/docs/sep92/hector1/paper_html/node2.html)

K:ev(1/k*matrix([k_x,0,0,0,k_z,k_y],[0,k_y,0,k_z,0,k_x],[0,0,k_z,k_y,k_x,0]),k_x=0,k_y=0);
k:ev(sqrt(k_x^2+k_y^2+k_z^2),k_x=0,k_y=0);

\begin{pmatrix} 0&0&0&0&{{\mathit{k}_{z}}\over{k}}&0\\ 0&0&0&{{\mathit{k}_{z}}\over{k}}&0&0\\ 0&0&{{\mathit{k}_{z}}\over{k}}&0&0&0 \end{pmatrix}

\left| \mathit{k}_{z}\right|

M:ratsimp(k^2/%rho*K.subst(C_33*(1+k_33^2),C_33,C).transpose(K));
ratsimp(eigenvalues(M));
ratsimp(eigenvectors(M));

\begin{pmatrix} {{C_{44}\,\mathit{k}_{z}^4}\over{\rho\,k^2}}&0&0\\ 0&{{C_{44}\,\mathit{k}_{z}^4}\over{\rho\,k^2}}&0\\ 0&0&{{\left(C_{33}\,k_{33}^2+C_{33}\right)\,\mathit{k}_{z}^4}\over{\rho\,k^2}} \end{pmatrix}

\left[ \left[ {{\left(C_{33}\,k_{33}^2+C_{33}\right)\,\mathit{k}_{z}^4}\over{\rho\,k^2}} , {{C_{44}\,\mathit{k}_{z}^4}\over{\rho\,k^2}} \right] , \left[ 1 , 2 \right] \right]

\left[ \left[ \left[ {{\left(C_{33}\,k_{33}^2+C_{33}\right)\,\mathit{k}_{z}^4}\over{\rho\,k^2}} , {{C_{44}\,\mathit{k}_{z}^4}\over{\rho\,k^2}} \right] , \left[ 1 , 2 \right] \right] , \left[ \left[ \left[ 0 , 0 , 1 \right] \right] , \left[ \left[ 1 , 0 , 0 \right] , \left[ 0 , 1 , 0 \right] \right] \right] \right]

The BAW vibration frequency, f_m, can be used to determine the thickness of the AlN plate, g_0

/* Without these, ratsimp and limit further down keep stopping to ask
   whether %epsilon_33 is positive, which makes the notebook impossible to
   batch and forces the reader to answer by hand.  0.2.4 declared
   only the four symbols on the first line. */
/* NB: %omega, %xi and %omega_m are deliberately NOT declared constant --
   limit(...,%omega,0) below needs %omega to be a variable.  The assume()
   calls are what silence the sign questions. */
declare([%rho,C_33,g_0,k_33,%epsilon_33,A,V],[constant,real,scalar])$
assume(%rho>0)$
assume(C_33>0)$
assume(g_0>0)$
assume(k_33>0)$
assume(%epsilon_33>0)$
assume(A>0)$
assume(%omega>0)$
assume(%xi>0)$
assume(%omega_m>0)$
assume(Q_m>0)$
k_z:2*%pi/(2*g_0);
%omega_m:sqrt(ev((k_z^2*C_33*(1+k_33^2))/(%rho)));
f_m:1/(2*%pi)*%omega_m;

{{\pi}\over{g_{0}}}

{{\sqrt{C_{33}}\,\pi\,\sqrt{k_{33}^2+1}}\over{\sqrt{\rho}\,g_{0}}}

{{\sqrt{C_{33}}\,\sqrt{k_{33}^2+1}}\over{2\,\sqrt{\rho}\,g_{0}}}

ratsimp(sqrt(%rho/(C_33*(1+k_33^2)))*%omega_m*g_0/2);

{{\pi}\over{2}}

m is the effective mass. A is the electrode area.

declare([A],[constant,real,scalar])$
/* Effective (modal) mass.  Two changes from 0.2.4:
     - integrate over the WHOLE plate, -g_0/2 .. g_0/2.  0.2.4 integrated
       over 0 .. g_0/2 only, i.e. half of it.
     - evaluate the mode shape at the fundamental thickness resonance,
       %xi = %pi/g_0, BEFORE integrating.  With %xi left symbolic Maxima
       stops and asks about floor(%xi*g_0/(4*%pi)) + ceiling(...) and never
       gets past it, which is one reason this notebook could not be batched.
   The answer is the familiar m = %rho*A*g_0/2, half the physical mass. */
m:A*%rho*integrate(ev(u_z(z)^2,%xi=%pi/g_0),z,-g_0/2,g_0/2)/A_1^2;

{{\rho\,A\,g_{0}\,\left(\epsilon_{33}\,C_{33}\,\xi\,g_{0}\,\cos \left({{\xi\,g_{0}}\over{2}}\right)\,k_{33}^2-2\,e_{33}^2\,\sin \left({{\xi\,g_{0}}\over{2}}\right)+\epsilon_{33}\,C_{33}\,\xi\,g_{0}\,\cos \left({{\xi\,g_{0}}\over{2}}\right)\right)^2}\over{2\,\left(-2\,e_{33}^2\,\sin \left({{g_{0}\,\xi}\over{2}}\right)+\epsilon_{33}\,C_{33}\,g_{0}\,k_{33}^2\,\xi\,\cos \left({{g_{0}\,\xi}\over{2}}\right)+\epsilon_{33}\,C_{33}\,g_{0}\,\xi\,\cos \left({{g_{0}\,\xi}\over{2}}\right)\right)^2}}

k_1_eff is the effective linear spring constant.

solve([%omega_m=sqrt(k_1_eff/m)],[k_1_eff])$
k_1_eff:ev(k_1_eff,%);

{{\pi^2\,A\,C_{33}\,k_{33}^2+\pi^2\,A\,C_{33}}\over{2\,g_{0}}}

Tf is the linear temperature coefficient of frequency.

/* Placeholder: nothing in f_m depends on Temp yet, so this evaluates to 0.
   Give the stiffness, density and thickness a temperature dependence to get
   a meaningful turnover curve. */
Tf:'diff(log(f_m),Temp);

{{d}\over{d\,\mathit{Temp}}}\,\log \left({{\sqrt{C_{33}}\,\sqrt{k_{33}^2+1}}\over{2\,\sqrt{\rho}\,g_{0}}}\right)

Charge, q, on the Electrodes

q:ratsimp(-D_z(g_0/2)*A);

-{{\left(\epsilon_{33}^2\,A\,C_{33}\,V\,k_{33}^2+\epsilon_{33}^2\,A\,C_{33}\,V\right)\,\xi\,\cos \left({{g_{0}\,\xi}\over{2}}\right)}\over{2\,e_{33}^2\,\sin \left({{g_{0}\,\xi}\over{2}}\right)+\left(-\epsilon_{33}\,C_{33}\,g_{0}\,k_{33}^2-\epsilon_{33}\,C_{33}\,g_{0}\right)\,\xi\,\cos \left({{g_{0}\,\xi}\over{2}}\right)}}

Motional Capacitance, C_m, Static Capacitance, C_0, and Transduction Factor, %eta

C_FBAR:ev(q/V,
    %xi=sqrt(%rho/(C_33*(1+e_33^2/(%epsilon_33*C_33))))*%omega,
    k_33=sqrt(e_33^2/(%epsilon_33*C_33)));
C_0:ratsimp(limit(C_FBAR,%omega,0));
ratsimp(C_FBAR-C_0)$
C_m:ratsimp(ev(k_33^2*C_0,
    k_33=sqrt(e_33^2/(%epsilon_33*C_33))));
/* Transduction factor: %eta = sqrt(C_m * k_1_eff), with C_m = k_33^2 C_0.
   0.2.4 carried an extra (1+k_33^2) factor here that QUARTZ.wxm did not; the
   two notebooks now agree. */
%eta:ratsimp(sqrt(C_m*k_1_eff));

-{{\sqrt{\rho}\,\left(\epsilon_{33}\,A\,V\,e_{33}^2+\epsilon_{33}^2\,A\,C_{33}\,V\right)\,\omega\,\cos \left({{\sqrt{\rho}\,g_{0}\,\omega}\over{2\,\sqrt{C_{33}}\,\sqrt{{{e_{33}^2}\over{\epsilon_{33}\,C_{33}}}+1}}}\right)}\over{\sqrt{C_{33}}\,V\,\sqrt{{{e_{33}^2}\over{\epsilon_{33}\,C_{33}}}+1}\,\left(2\,e_{33}^2\,\sin \left({{\sqrt{\rho}\,g_{0}\,\omega}\over{2\,\sqrt{C_{33}}\,\sqrt{{{e_{33}^2}\over{\epsilon_{33}\,C_{33}}}+1}}}\right)+{{\sqrt{\rho}\,\left(-e_{33}^2\,g_{0}-\epsilon_{33}\,C_{33}\,g_{0}\right)\,\omega\,\cos \left({{\sqrt{\rho}\,g_{0}\,\omega}\over{2\,\sqrt{C_{33}}\,\sqrt{{{e_{33}^2}\over{\epsilon_{33}\,C_{33}}}+1}}}\right)}\over{\sqrt{C_{33}}\,\sqrt{{{e_{33}^2}\over{\epsilon_{33}\,C_{33}}}+1}}}\right)}}

{{A\,e_{33}^2+\epsilon_{33}\,A\,C_{33}}\over{C_{33}\,g_{0}}}

{{A\,e_{33}^4+\epsilon_{33}\,A\,C_{33}\,e_{33}^2}\over{\epsilon_{33}\,C_{33}^2\,g_{0}}}

{{\sqrt{A\,e_{33}^4+\epsilon_{33}\,A\,C_{33}\,e_{33}^2}\,\sqrt{\pi^2\,A\,C_{33}\,k_{33}^2+\pi^2\,A\,C_{33}}}\over{\sqrt{2}\,\sqrt{\epsilon_{33}}\,C_{33}\,g_{0}}}

Material Parameters

Hexagonal Aluminum Nitride (Wurtzite, Crystal Class 6mm)

(See "Material Parameters of AlN and LiAlO2 Single Crystals", by A. V. Sotnikov, H. Schmidt, M. Weihnacht, E. P. Smirnova, T. Yu. Chemekova, Yu. N. Makarov, IEEE IFCS, pp. 935-938, 2009)

C_11:402.5;
C_12:135.6;
C_13:101;
C_33:387.6;
C_44:122.9;
C:1e9*matrix([C_11, C_12, C_13, 0, 0, 0],[C_12, C_11, C_13, 0, 0, 0],[C_13, C_13, C_33, 0, 0, 0],[0, 0, 0, C_44, 0, 0],[0, 0, 0, 0, C_44, 0],[0, 0, 0, 0, 0, (C_11-C_12)/2]);

402.5

135.6

101

387.6

122.9

\begin{pmatrix} 4.025 \times 10^{+11}&1.356 \times 10^{+11}&1.01 \times 10^{+11}&0&0&0\\ 1.356 \times 10^{+11}&4.025 \times 10^{+11}&1.01 \times 10^{+11}&0&0&0\\ 1.01 \times 10^{+11}&1.01 \times 10^{+11}&3.876 \times 10^{+11}&0&0&0\\ 0&0&0&1.229 \times 10^{+11}&0&0\\ 0&0&0&0&1.229 \times 10^{+11}&0\\ 0&0&0&0&0&1.3345 \times 10^{+11} \end{pmatrix}

e_33:1.34;
e_31:-0.6;
e_15:-0.32;
e:matrix([0,0,0,0,e_15,0],[0,0,0,e_15,0,0],[e_31,e_31,e_33,0,0,0]);

1.34

-0.6

-0.32

\begin{pmatrix} 0&0&0&0&-0.32&0\\ 0&0&0&-0.32&0&0\\ -0.6&-0.6&1.34&0&0&0 \end{pmatrix}

%epsilon_0:8.854e-12;
%epsilon_11:9;
%epsilon_33:9.5;
%epsilon:%epsilon_0*matrix([%epsilon_11,0,0],[0,%epsilon_11,0],[0,0,%epsilon_33]);

8.854 \times 10^{-12}

9

9.5

\begin{pmatrix} 7.968599999999999 \times 10^{-11}&0&0\\ 0&7.968599999999999 \times 10^{-11}&0\\ 0&0&8.4113 \times 10^{-11} \end{pmatrix}

%rho:3260;

3260