Skip to content
Snippets Groups Projects
Commit bdef493f authored by Sam Katiraee-Far's avatar Sam Katiraee-Far
Browse files

Added meeting notebook for today

parent 1725cb03
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Meeting Notes/Meeting 06_07_2021.ipynb 0 → 100644
+ 64
− 0
View file @ bdef493f
%% Cell type:markdown id: tags:
# Agenda
* Look at results from
* Wave equation
* Looks good to me
* 1D Homogeneous
* Looks good to me
* 1D space dependent
* Matrix of space independent case is already symmetric
* Do not get the same result when I apply the specific homogeneous case to the general space dependent method
* Matrix has same shape as homogeneous notebook, but different values
* Gives same modes
* But different frequencies
* Last night noticed wrong minus sign in 2nd derivative FD :(
* However, I don't think that effects this specific problem
* Boundary conditions
* is du/dx = 0 at boundaries correct?
* How to move forward
* Fix the minus sign
* Good to do, but I don't think it will solve the issue
* Have to think about what would solve it
* If this version works ->
* Apply to some space dependent coefficient cases
* Which ones would we want?
* How to check correctness of results?
* There also is a 1D anisotropic case in Tiersten paper
* Study waves in 1D, but allow for anisotropic coefficients
* Time dependent?
* 2D?
Simulations/1D_FD_SpaceDep.ipynb
+ 21
− 17
View file @ bdef493f
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pylab as plt
from time import time
import sys
import scipy.integrate as sp
```
%% Cell type:markdown id: tags:
We start with piezo-electric coupled equations
\begin{equation}
T = c^E \frac{\partial u}{\partial x} - e E \qquad (1) \\
D = e \frac{\partial u}{\partial x} + \epsilon^S E \qquad (2)
\end{equation}
where the physical quantities denoted by the symbols is given in the table below
| Symbol | Quanitity |
| ----------------|-------------|
| $T \thinspace \mathrm{(Pa)}$ | Stress |
| $E \thinspace \mathrm{(V/m)}$ | Electric field |
| $c^E \thinspace \mathrm{(Pa)}$ | Stiffness coefficient |
| $u \thinspace \mathrm{(m)}$ | Mechanical displacement |
| $e \thinspace \mathrm{(C/m^2)}$ | Piezoelectric coefficient |
| $D \thinspace \mathrm{(C/m^2)}$ | Electric displacement |
| $\epsilon^S \thinspace \mathrm{(C/Vm)}$ | Permitivitty under constant strain $S$ |
We can use the equations
$$
\frac{\partial D}{\partial x} = 0
$$
and the wave equation for acoustic media
$$
\frac{\partial T}{\partial x} = \rho \frac{\partial^2 u}{\partial t^2 }
$$
combined with equation (1) and (2) to obtain
$$
\rho \frac{\partial^2 u}{\partial t^2 }= \frac{\partial c^E }{\partial x} \frac{\partial u}{\partial x} + c^E \frac{\partial^2 u}{\partial x^2}
+ \frac{\partial e}{\partial x} \frac{\partial \phi}{\partial x} + e \frac{\partial^2 \phi}{\partial x^2} \qquad (3)\\
\frac{\partial e}{\partial x} \frac{\partial u}{\partial x} + e \frac{\partial^2 u}{\partial x^2} - \left(\frac{\partial \epsilon^S}{\partial x} \frac{\partial \phi}{\partial x} + \epsilon^S \frac{\partial^2 \phi}{\partial x^2} \right)= 0 \qquad (4)\\
$$
where it was used that the electric field is minus the gradient of the potential $\phi$
and that all quantities, including constants and coefficients, are dependent on $x$.
We can now rewrite equations (3) and (4) using the finite difference approximation, where we use the central difference formula. For the first derivative this is given by
$$
\frac{\partial f_i}{\partial x} \approx \frac{f_{i+1}-f_{i-1}}{2h},
$$
and for the second derivative it is given as
$$
-\frac{\partial^2 f_i}{\partial x^2} \approx \frac{-f_{i+1} + 2f_{i} -f_{i-1}}{h^2}
\frac{\partial^2 f_i}{\partial x^2} \approx \frac{f_{i+1} - 2f_{i} +f_{i-1}}{h^2}
$$
where $h$ is the discretization size and $f_i$ is the value of $f(x)$ at the $i^{th}$ point.
Using this we obtain the finite difference equations of equations (3) and (4):
$$
\rho_i \frac{\partial^2 u_i}{\partial t^2 }= c^E_{i,x} \frac{u_{i+1}-u_{i-1}}{2h} + c^E_i \frac{-u_{i+1}+2u_i-u_{i-1}}{h^2}
+ e_{i,x} \frac{\phi_{i+1}-\phi_{i-1}}{2h} + e_i \frac{-\phi_{i+1}+2\phi_{i}-\phi_{i-1}}{h^2} \qquad (5)\\
e_{i,x} \frac{u_{i+1}-u_{i-1}}{2h} + e_i \frac{-u_{i+1}+2u_i-u_{i-1}}{h^2} - \left(\epsilon_{i,x} \frac{\phi_{i+1}-\phi_{i-1}}{2h} + \epsilon_i^S \frac{-\phi_{i+1}+2\phi_{i}-\phi_{i-1}}{h^2} \right)= 0 \qquad (6)\\
\rho_i \frac{\partial^2 u_i}{\partial t^2 }= c^E_{i,x} \frac{u_{i+1}-u_{i-1}}{2h} + c^E_i \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}
+ e_{i,x} \frac{\phi_{i+1}-\phi_{i-1}}{2h} + e_i \frac{\phi_{i+1}-2\phi_{i}+\phi_{i-1}}{h^2} \qquad (5)\\
e_{i,x} \frac{u_{i+1}-u_{i-1}}{2h} + e_i \frac{u_{i+1}-2u_i+u_{i-1}}{h^2} - \epsilon_{i,x} \frac{\phi_{i+1}-\phi_{i-1}}{2h} - \epsilon_i^S \frac{\phi_{i+1}-2\phi_{i}+\phi_{i-1}}{h^2} = 0 \qquad (6)\\
$$
where we denote the derivative with respect to $x$ of the coefficients $c^E, e, \epsilon^S$ by placing it after the comma. We can rewrite equations (5) and (6) in a form where all the discretized variables only appear no more than once on each side of the equation:
$$
\rho_i \frac{\partial^2 u_i}{\partial t^2 } = u_{i-1} \left(-\frac{c^E_{i,x}}{2h} - \frac{c^E_i}{h^2}\right) + u_i \left( \frac{2c^E_i}{h^2} \right) + u_{i+1} \left(\frac{c^E_{i,x}}{2h} - \frac{c^E_i}{h^2}\right) + \phi_{i-1} \left(- \frac{e_{i,x}}{2h} - \frac{e_i}{h^2} \right) + \phi_i \left(\frac{2e_i}{h^2}\right) + \phi_{i+1} \left( \frac{e_{i,x}}{2h} - \frac{e_i}{h^2} \right) \quad(7)\\
u_{i-1} \left( -\frac{e_{i,x}}{2h} - \frac{e_i}{h^2} \right) + u_i \left(\frac{2e_i}{h^2} \right) + u_{i+1} \left( \frac{e_{i,x}}{2h} - \frac{e_i}{h^2} \right) + \phi_{i-1} \left(+\frac{\epsilon^S_{i,x}}{2h} + \frac{\epsilon^S_i}{h^2} \right) + \phi_i \left(-\frac{2\epsilon^S_i}{h^2}\right) + \phi_{i+1} \left( -\frac{\epsilon^S_{i,x}}{2h} + \frac{\epsilon^S_i}{h^2} \right) = 0 \qquad(8)
\rho_i \frac{\partial^2 u_i}{\partial t^2 } = u_{i-1} \left(-\frac{c^E_{i,x}}{2h} + \frac{c^E_i}{h^2}\right) + u_i \left( -\frac{2c^E_i}{h^2} \right) + u_{i+1} \left(\frac{c^E_{i,x}}{2h} +\frac{c^E_i}{h^2}\right) + \phi_{i-1} \left(- \frac{e_{i,x}}{2h} + \frac{e_i}{h^2} \right) + \phi_i \left(-\frac{2e_i}{h^2}\right) + \phi_{i+1} \left( \frac{e_{i,x}}{2h} + \frac{e_i}{h^2} \right) \quad(7)\\
u_{i-1} \left( -\frac{e_{i,x}}{2h} + \frac{e_i}{h^2} \right) + u_i \left(-\frac{2e_i}{h^2} \right) + u_{i+1} \left( \frac{e_{i,x}}{2h} + \frac{e_i}{h^2} \right) + \phi_{i-1} \left(+\frac{\epsilon^S_{i,x}}{2h} - \frac{\epsilon^S_i}{h^2} \right) + \phi_i \left(\frac{2\epsilon^S_i}{h^2}\right) + \phi_{i+1} \left( -\frac{\epsilon^S_{i,x}}{2h} - \frac{\epsilon^S_i}{h^2} \right) = 0 \qquad(8)
$$
We use a free boundary condition for the mechanical displacement, which means that its derivative is zero there. Numerically this means that:
$$
u_{N-1} - u_{N-2}=0 \quad \rightarrow \quad u_{N-1} = u_{N-2}\\
u_1 - u_{0} = 0 \quad \rightarrow \quad u_0 = u_1 \\
$$
For the potential we take that its value at the boundary is 0
$$
\phi_0 = 0 \\
\phi_{N-1} = 0 \\
$$
%% Cell type:markdown id: tags:
So the equations we obtain for the points near the boundary $i = 1,2, N-3, N-2$ are
$i = 1$:
$$
\frac{\partial^2 u_1}{\partial t^2 } = u_{1} \left(-\frac{c^E_{1,x}}{2h} + \frac{c^E_1}{h^2} \right) + u_{2} \left(\frac{c^E_{1,x}}{2h} - \frac{c^E_1}{h^2}\right) + \phi_1 \left(\frac{2e_1}{h^2}\right) + \phi_{2} \left( \frac{e_{1,x}}{2h} - \frac{e_1}{h^2} \right) \\
u_{1} \left( -\frac{e_{1,x}}{2h} - \frac{e_1}{h^2} + \frac{2e_1}{h^2} \right) + u_{2} \left( \frac{e_{1,x}}{2h} - \frac{e_1}{h^2} \right) + \phi_1 \left(-\frac{2\epsilon^S_1}{h^2}\right) + \phi_{2} \left( -\frac{\epsilon^S_{1,x}}{2h} + \frac{\epsilon^S_1}{h^2} \right) = 0
\frac{\partial^2 u_1}{\partial t^2 } = u_{1} \left(-\frac{c^E_{1,x}}{2h} - \frac{c^E_1}{h^2} \right) + u_{2} \left(\frac{c^E_{1,x}}{2h} + \frac{c^E_1}{h^2}\right) + \phi_1 \left(-\frac{2e_1}{h^2}\right) + \phi_{2} \left( \frac{e_{1,x}}{2h} + \frac{e_1}{h^2} \right) \\
u_{1} \left( -\frac{e_{1,x}}{2h} - \frac{e_1}{h^2} \right) + u_{2} \left( \frac{e_{1,x}}{2h} + \frac{e_1}{h^2} \right) + \phi_1 \left(\frac{2\epsilon^S_1}{h^2}\right) + \phi_{2} \left( -\frac{\epsilon^S_{1,x}}{2h} - \frac{\epsilon^S_1}{h^2} \right) = 0
$$
$i=2$:
$$
\frac{\partial^2 u_2}{\partial t^2 } = u_{1} \left(-\frac{c^E_{2,x}}{2h} - \frac{c^E_2}{h^2}\right) + u_2 \left( \frac{2c^E_2}{h^2} \right) + u_{3} \left(\frac{c^E_{2,x}}{2h} - \frac{c^E_2}{h^2}\right) + \phi_{1} \left(- \frac{e_{2,x}}{2h} - \frac{e_2}{h^2} \right) + \phi_2 \left(\frac{2e_2}{h^2}\right) + \phi_{3} \left( \frac{e_{2,x}}{2h} - \frac{e_2}{h^2} \right) \\
u_{1} \left( -\frac{e_{2,x}}{2h} - \frac{e_2}{h^2} \right) + u_2 \left(\frac{2e_2}{h^2} \right) + u_{3} \left( \frac{e_{2,x}}{2h} - \frac{e_2}{h^2} \right) + \phi_{1} \left(+\frac{\epsilon^S_{2,x}}{2h} + \frac{\epsilon^S_2}{h^2} \right) + \phi_2 \left(-\frac{2\epsilon^S_2}{h^2}\right) + \phi_{3} \left( -\frac{\epsilon^S_{2,x}}{2h} + \frac{\epsilon^S_2}{h^2} \right) = 0
\frac{\partial^2 u_2}{\partial t^2 } = u_{1} \left(-\frac{c^E_{2,x}}{2h} + \frac{c^E_2}{h^2}\right) + u_2 \left( -\frac{2c^E_2}{h^2} \right) + u_{3} \left(\frac{c^E_{2,x}}{2h} + \frac{c^E_2}{h^2}\right) + \phi_{1} \left(- \frac{e_{2,x}}{2h} + \frac{e_2}{h^2} \right) + \phi_2 \left(-\frac{2e_2}{h^2}\right) + \phi_{3} \left( \frac{e_{2,x}}{2h} + \frac{e_2}{h^2} \right) \\
u_{1} \left( -\frac{e_{2,x}}{2h} + \frac{e_2}{h^2} \right) + u_2 \left(-\frac{2e_2}{h^2} \right) + u_{3} \left( \frac{e_{2,x}}{2h} + \frac{e_2}{h^2} \right) + \phi_{1} \left(\frac{\epsilon^S_{2,x}}{2h} - \frac{\epsilon^S_2}{h^2} \right) + \phi_2 \left(\frac{2\epsilon^S_2}{h^2}\right) + \phi_{3} \left( -\frac{\epsilon^S_{2,x}}{2h} - \frac{\epsilon^S_2}{h^2} \right) = 0
$$
$i=N-3$
$$
\frac{\partial^2 u_{N-3}}{\partial t^2 } = u_{N-4} \left(-\frac{c^E_{N-3,x}}{2h} - \frac{c^E_{N-3}}{h^2}\right) + u_i \left( \frac{2c^E_{N-3}}{h^2} \right) + u_{N-2} \left(\frac{c^E_{N-3,x}}{2h} - \frac{c^E_{N-3}}{h^2}\right) + \phi_{N-4} \left(- \frac{e_{N-4,x}}{2h} - \frac{e_{N-3}}{h^2} \right) + \phi_{N-3} \left(\frac{2e_{N-3}}{h^2}\right) + \phi_{N-2} \left( \frac{e_{N-3,x}}{2h} - \frac{e_{N-3}}{h^2} \right) \\
u_{N-4} \left( -\frac{e_{N-3,x}}{2h} - \frac{e_N-3}{h^2} \right) + u_{N-3} \left(\frac{2e_{N-3}}{h^2} \right) + u_{N-2} \left( \frac{e_{N-3,x}}{2h} - \frac{e_{N-3}}{h^2} \right) + \phi_{N-4} \left(+\frac{\epsilon^S_{N-3,x}}{2h} + \frac{\epsilon^S_{N-3}}{h^2} \right) + \phi_{N-3} \left(-\frac{2\epsilon^S_{N-3}}{h^2}\right) + \phi_{N-2} \left( -\frac{\epsilon^S_{{N-3},x}}{2h} + \frac{\epsilon^S_{N-3}}{h^2} \right) = 0
\frac{\partial^2 u_{N-3}}{\partial t^2 } = u_{N-4} \left(-\frac{c^E_{N-3,x}}{2h} + \frac{c^E_{N-3}}{h^2}\right) + u_{N-3} \left( - \frac{2c^E_{N-3}}{h^2} \right) + u_{N-2} \left(\frac{c^E_{N-3,x}}{2h} + \frac{c^E_{N-3}}{h^2}\right) + \phi_{N-4} \left(- \frac{e_{N-4,x}}{2h} + \frac{e_{N-3}}{h^2} \right) + \phi_{N-3} \left(-\frac{2e_{N-3}}{h^2}\right) + \phi_{N-2} \left( \frac{e_{N-3,x}}{2h} + \frac{e_{N-3}}{h^2} \right) \\
u_{N-4} \left( -\frac{e_{N-3,x}}{2h} + \frac{e_N-3}{h^2} \right) + u_{N-3} \left(-\frac{2e_{N-3}}{h^2} \right) + u_{N-2} \left( \frac{e_{N-3,x}}{2h} + \frac{e_{N-3}}{h^2} \right) + \phi_{N-4} \left(\frac{\epsilon^S_{N-3,x}}{2h} - \frac{\epsilon^S_{N-3}}{h^2} \right) + \phi_{N-3} \left(\frac{2\epsilon^S_{N-3}}{h^2}\right) + \phi_{N-2} \left( -\frac{\epsilon^S_{{N-3},x}}{2h} - \frac{\epsilon^S_{N-3}}{h^2} \right) = 0
$$
$i=N-2$
$$
\frac{\partial^2 u_{N-2}}{\partial t^2 } = u_{N-3} \left(-\frac{c^E_{N-2,x}}{2h} - \frac{c^E_N-2}{h^2}\right) + u_{N-2} \left( \frac{2c^E_i}{h^2} + \frac{c^E_{N-2,x}}{2h} - \frac{c^E_{N-2}}{h^2}\right) + \phi_{N-3} \left(- \frac{e_{N-2,x}}{2h} - \frac{e_N-2}{h^2} \right) + \phi_{N-2} \left(\frac{2e_i}{h^2}\right) \\
\frac{\partial^2 u_{N-2}}{\partial t^2 } = u_{N-3} \left(-\frac{c^E_{N-2,x}}{2h} - \frac{c^E_{N-2}}{h^2}\right) + u_{N-2} \left( \frac{2c^E_i}{h^2} + \frac{c^E_{N-2,x}}{2h} - \frac{c^E_{N-2}}{h^2}\right) + \phi_{N-3} \left(- \frac{e_{N-2,x}}{2h} - \frac{e_N-2}{h^2} \right) + \phi_{N-2} \left(\frac{2e_i}{h^2}\right) \\
u_{N-3} \left( -\frac{e_{N-2,x}}{2h} - \frac{e_{N-2}}{h^2} \right) + u_{N-2} \left(\frac{e_{N-2}}{h^2} + \frac{e_{N-2,x}}{2h} \right) + \phi_{N-3} \left(\frac{\epsilon^S_{N-2,x}}{2h} + \frac{\epsilon^S_{N-2}}{h^2} \right) + \phi_{N-2} \left(-\frac{2\epsilon^S_{N-2}}{h^2}\right) = 0
$$
%% Cell type:markdown id: tags:
We can simplify assuming that we can seperate the variables as $u(x,t) = u(x)e^{-i \omega t}$. We then obtain
$$
-\rho_i \omega^2 u_i = u_{i-1} \left(-\frac{c^E_{i,x}}{2h} - \frac{c^E_i}{h^2}\right) + u_i \left( \frac{2c^E_i}{h^2} \right) + u_{i+1} \left(\frac{c^E_{i,x}}{2h} - \frac{c^E_i}{h^2}\right) + \phi_{i-1} \left(- \frac{e_{i,x}}{2h} - \frac{e_i}{h^2} \right) + \phi_i \left(\frac{2e_i}{h^2}\right) + \phi_{i+1} \left( \frac{e_{i,x}}{2h} - \frac{e_i}{h^2} \right) \quad \\
u_{i-1} \left( -\frac{e_{i,x}}{2h} - \frac{e_i}{h^2} \right) + u_i \left(\frac{2e_i}{h^2} \right) + u_{i+1} \left( \frac{e_{i,x}}{2h} - \frac{e_i}{h^2} \right) + \phi_{i-1} \left(+\frac{\epsilon^S_{i,x}}{2h} + \frac{\epsilon^S_i}{h^2} \right) + \phi_i \left(-\frac{2\epsilon^S_i}{h^2}\right) + \phi_{i+1} \left( -\frac{\epsilon^S_{i,x}}{2h} + \frac{\epsilon^S_i}{h^2} \right) = 0
$$
%% Cell type:markdown id: tags:
These equations can be rewritten in matrix form
$$
- \omega^2 \pmb{\rho} \circ \mathbf u = A \mathbf{u} + B \pmb{\phi} \quad (9) \\
C \mathbf{u} + D\pmb{\phi} = \mathbf{0} \qquad (10)\\
$$
%% Cell type:markdown id: tags:
where the matrices $A, B, C, D$ are given as:
$$
A = \begin{pmatrix}
-\frac{c^E_{1,x}}{2h} + \frac{c^E_1}{h^2} & \frac{c^E_{1,x}}{2h} - \frac{c^E_1}{h^2} & 0 & 0&\dots & 0 & 0 & 0&0\\
-\frac{c^E_{2,x}}{2h} - \frac{c^E_2}{h^2} & \frac{2c^E_2}{h^2} & \frac{c^E_{2,x}}{2h} - \frac{c^E_2}{h^2} & 0 & \dots & 0 & 0 & 0&0\\
&&&& \ddots &&&&\\
0&0&0&0& \dots &0 &-\frac{c^E_{N-3,x}}{2h} - \frac{c^E_{N-3}}{h^2} & \frac{2c^E_{N-3}}{h^2} & \frac{c^E_{N-3,x}}{2h} - \frac{c^E_{N-3}}{h^2}\\
0&0&0&0&\dots&0&0& -\frac{c^E_{N-2,x}}{2h} - \frac{c^E_{N-2}}{h^2} & \frac{c^E_{N-2}}{h^2} + \frac{c^E_{N-2,x}}{2h}\\
\end{pmatrix}
$$
$$
B = \begin{pmatrix}
\frac{2e_1}{h^2} & \frac{e_{1,x}}{2h} - \frac{e_1}{h^2} & 0 & 0 &\dots & 0 & 0 & 0 & 0 \\
- \frac{e_{2,x}}{2h} - \frac{e_2}{h^2} & \frac{2e_2}{h^2} & \frac{e_{2,x}}{2h} - \frac{e_2}{h^2} & 0 &\dots & 0 & 0 & 0 & 0 \\
&&&& \ddots &&&&\\
0 & 0 & 0 & 0 & \dots & 0 & - \frac{e_{N-3,x}}{2h} - \frac{e_{N-3}}{h^2} & \frac{2e_{N-3}}{h^2} & \frac{e_{N-3,x}}{2h} - \frac{e_{N-3}}{h^2} \\
0&0&0&0&\dots&0&0& - \frac{e_{N-2,x}}{2h} - \frac{e_N-2}{h^2} & \frac{2e_i}{h^2}\\
\end{pmatrix}
$$
$$
C = \begin{pmatrix}
-\frac{e_{1,x}}{2h} + \frac{e_1}{h^2} & \frac{e_{1,x}}{2h} - \frac{e_1}{h^2} &0&0&\dots&0&0&0&0\\
-\frac{e_{2,x}}{2h} - \frac{e_2}{h^2} & \frac{2e_2}{h^2} & \frac{e_{2,x}}{2h} - \frac{e_2}{h^2} &0 & \dots &0&0&0&0 \\
&&&& \ddots &&&&\\
0&0&0&0&\dots&0& -\frac{e_{N-3,x}}{2h} - \frac{e_{N-3}}{h^2} & \frac{2e_{N-3}}{h^2} & \frac{e_{N-3,x}}{2h} - \frac{e_{N-3}}{h^2}\\
0&0&0&0&\dots&0&0& -\frac{e_{N-2,x}}{2h} - \frac{e_{N-2}}{h^2} & \frac{e_{N-2}}{h^2} + \frac{e_{N-2,x}}{2h}\\
\end{pmatrix}
$$
$$
D = \begin{pmatrix}
-\frac{2\epsilon^S_1}{h^2} & -\frac{\epsilon^S_{1,x}}{2h} + \frac{\epsilon^S_1}{h^2} & 0 & 0 & \dots & 0 & 0 & 0 &0 \\
\frac{\epsilon^S_{2,x}}{2h} + \frac{\epsilon^S_2}{h^2} & -\frac{2\epsilon^S_2}{h^2} & -\frac{\epsilon^S_{2,x}}{2h} + \frac{\epsilon^S_2}{h^2} & 0 &\dots &0&0&0&0 \\
&&&& \ddots &&&&\\
0&0&0&0&\dots&0& \frac{\epsilon^S_{N-3,x}}{2h} + \frac{\epsilon^S_{N-3}}{h^2} & -\frac{2\epsilon^S_{N-3}}{h^2}& -\frac{\epsilon^S_{{N-3},x}}{2h} + \frac{\epsilon^S_{N-3}}{h^2}\\
0&0&0&0&\dots&0&0& \frac{\epsilon^S_{N-2,x}}{2h} + \frac{\epsilon^S_{N-2}}{h^2} & -\frac{2\epsilon^S_{N-2}}{h^2} \\
\end{pmatrix}
$$
This is the generalized eigenvalue problem we want to solve. We can then extract $\pmb{\phi}$ from equation (10):
$$
\pmb{\phi} = -D^{-1} C \mathbf{u}
$$
and substitute it into equation (9), to obtain:
$$
\omega^2 \mathbf{u} = (BD^{-1}C- A) \mathbf{u} \circ \pmb{\rho}^{-1}
$$
%% Cell type:markdown id: tags:
# Code
%% Cell type:code id: tags:
``` python
#compute eigenfrequencies and modes
def find_solution(interval, N, eps_r, rho, eps, eps_p, e, e_p, c, c_p):
'''
Input
----------------------
interval: tuple
beginning and end coordinates
N: integer
number of points in discretization
eps_r: float
relative permitivitty
rho: callable
Density. Should take parameter x
eps: callable
permitivitty dependent on space. Should take parameter x and relative permitivitty eps_r
eps_p: callable
derivative of permitivitty. Should take parameter x and relative permitivitty eps_r
e: callable
piezo electric coefficient. Should take paremeter x.
e_p: callable
derivative of piezo electric coefficient. Should take parameter x.
c: callable
stiffness coefficient. Should take parameter x.
c_p: callable
Derivative of stiffness coefficient. Should take parameter x.
Output
--------------------
Discretization,E igenfrequencies, Eigenmodes
'''
#vectorize the input functions
eps = np.vectorize(eps)
eps_p = np.vectorize(eps_p)
e = np.vectorize(e)
e_p = np.vectorize(e_p)
c = np.vectorize(c)
c_p = np.vectorize(c_p)
rho = np.vectorize(rho)
#make x array
x0, x1 = interval
x = np.linspace(x0,x1,N)
dx = x[1] - x[0]
#make the arrays with coefficients
rho_arr = rho(x[1:-1])
eps_arr = eps(eps_r,x[1:-1])
eps_p_arr = eps_p(eps_r,x[1:-1])
e_arr = e(x[1:-1])
e_p_arr = e_p(x[1:-1])
c_arr = c(x[1:-1])
c_p_arr = c_p(x[1:-1])
#Make the matrices
#matrix A
A0 = 2 * c_arr / dx**2 #main diagonal
A0[0] = c_arr[0] / dx**2 - c_p_arr[0] / (2*dx) #first point B.C.
A0[-1] = c_arr[-1] / dx**2 - c_p_arr[-1] / (2*dx) #first point B.C. #last point B.C.
AR = c_p_arr[:-1] / (2*dx) - c_arr[:-1] / (dx**2) #Right diagonal
AL = -c_p_arr[1:] / (2*dx) - c_arr[1:] / (dx**2) #left diagonal
A = np.diag(A0) + np.diag(AR, k=1) + np.diag(AL, k=-1)
#matrix C (similar to A)
C0 = 2 * e_arr / dx**2 #main diagonal
C0[0] = e_arr[0] / dx**2 - e_p_arr[0] / (2*dx) #first point B.C.
C0[-1] = e_arr[-1] / dx**2 - e_p_arr[-1] / (2*dx) #last point B.C.
CR = e_p_arr[:-1] / (2*dx) - e_arr[:-1] / (dx**2) #Right diagonal
CL = -e_p_arr[1:] / (2*dx) - e_arr[1:] / (dx**2) #left diagonal
C = np.diag(C0) + np.diag(CR, k=1) + np.diag(CL, k=-1)
#matrix B
B0 = 2 * e_arr / dx**2 #main diagonal
BR = e_p_arr[:-1] / (2*dx) - e_arr[:-1] / (dx**2) #Right diagonal
BL = -e_p_arr[1:] / (2*dx) - e_arr[1:] / (dx**2) #left diagonal
B = np.diag(B0) + np.diag(BR, k=1) + np.diag(BL, k=-1)
#matrix D (similar to B)
D0 = -2 * e_arr / dx**2 #main diagonal
DR = -e_p_arr[:-1] / (2*dx) + e_arr[:-1] / (dx**2) #Right diagonal
DL = e_p_arr[1:] / (2*dx) + e_arr[1:] / (dx**2) #left diagonal
D = np.diag(D0) + np.diag(DR, k=1) + np.diag(DL, k=-1)
Dinv = np.linalg.inv(D)
#Find the matrix of the eigenvalue problem
M = 1/rho_arr * (A - np.dot(B,np.dot(Dinv,C)))
# print("A")
# print(np.round(A, decimals = 3))
# print("B")
# print(np.round(B, decimals = 3))
# print("C")
# print(np.round(C, decimals = 3))
# print("D")
# print(np.round(D, decimals = 3))
np.set_printoptions(precision=1, threshold=sys.maxsize)
print("M")
print(np.round(M, decimals = 0))
#find the squared eigenvalues and eigenvectors
eigvals, eigvectors = np.linalg.eigh(M)
#sort the eigenfrequencies
eigfreq_sorted = np.sort(eigvals)
#sort the eigenmodes
eigvectors_sorted = eigvectors[:,np.argsort(eigvals)]
#duplicate first and last rows to take into account the zero neumann B.C.
v_first = np.array([eigvectors_sorted[0]])
v_last = np.array([eigvectors_sorted[-1]])
eigvectors_sorted = np.concatenate((v_first,eigvectors_sorted,v_last))
return x, eigfreq_sorted, eigvectors_sorted
```
%% Cell type:markdown id: tags:
# Test method by comparing to homogeneous case
Using parameters for AlN. Setting all the derivatives to zero.
range: 0 - 200$\mu$m (Jaspers HBAR)
from thesis <font color = "blue">Bulk Acoustic Wave Resonators and their Application to Microwave Devices<font color = "black">:
* $e$ = 1.5 C/m^2
* $c$ = 395 N/m^2
* $\rho$ = 3260 kg/m^3
* $\epsilon_r$ = 10.5
All derivatives are zero as we study the homogeneous case
%% Cell type:code id: tags:
``` python
#WHAT DO WE NEED?
#Functions which describe the coefficients dependent on space
#Will take constant value to check first
def permitivitty(eps_r, x):
eps0 = 8.85419e-12 #[F/m] (rounded to 6 singificant digits)
eps = eps_r * eps0
return eps
def permitivitty_prime(eps_r, x):
eps0 = 8.85419e-12 #[F/m] (rounded to 6 singificant digits)
return 0
def piezoCoefficient(x):
e = 1.5 #C/m^2
return e
def piezoCoefficient_prime(x):
return 0
def stiffness(x):
return 395 #[Pa]
def stiffness_prime(x):
return 0
def density(x):
return 3260 #[kg/m^3]
```
%% Cell type:code id: tags:
``` python
#test the function
x0 = 0
x1 = 200e-6
N = 10
eps_r = 10.5
x, eigvals, eigmodes = find_solution([x0,x1],N,eps_r,density,permitivitty, permitivitty_prime, \
piezoCoefficient, piezoCoefficient_prime, stiffness, stiffness_prime)
#analytic frequency calculation
e = piezoCoefficient(1)
eps = permitivitty(10.5,1)
c = stiffness(1)
rho = density(1)
k = np.arange(1,N+1,1)
omega_analytic = k * np.pi / (x1-x0) * np.sqrt((c + e/eps)/rho)
eigfreqs_sorted = np.sqrt((np.sort(eigvals)))
plt.plot(eigfreqs_sorted)
#plt.plot(omega_analytic)
#plot figure
plt.figure(figsize=(14,4))
plt.plot(x, eigmodes[:,3])
plt.xlabel("x (m)")
plt.ylabel("u (m)")
```
%% Output
M
[[ 2.5e+08 -2.5e+08 0.0e+00 -0.0e+00 -0.0e+00 0.0e+00 0.0e+00 0.0e+00]
[-2.5e+08 4.9e+08 -2.5e+08 0.0e+00 -0.0e+00 0.0e+00 -0.0e+00 0.0e+00]
[-0.0e+00 -2.5e+08 4.9e+08 -2.5e+08 0.0e+00 -0.0e+00 0.0e+00 -0.0e+00]
[ 0.0e+00 -0.0e+00 -2.5e+08 4.9e+08 -2.5e+08 0.0e+00 -0.0e+00 0.0e+00]
[-0.0e+00 0.0e+00 -0.0e+00 -2.5e+08 4.9e+08 -2.5e+08 0.0e+00 -0.0e+00]
[ 0.0e+00 -0.0e+00 -0.0e+00 0.0e+00 -2.5e+08 4.9e+08 -2.5e+08 0.0e+00]
[ 0.0e+00 0.0e+00 0.0e+00 -0.0e+00 0.0e+00 -2.5e+08 4.9e+08 -2.5e+08]
[ 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 -2.5e+08 2.5e+08]]
/Users/samkatiraee-far/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:20: RuntimeWarning: invalid value encountered in sqrt
Text(0,0.5,'u (m)')
%% Cell type:code id: tags:
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment