commit 54f4f0161af99a7e057e6385981e5010b4da58c2
parent 466fd9b5e85c1daae1be12b34431d4dc56e1d816
Author: jatinchowdhury18 <[email protected]>
Date: Wed, 6 Feb 2019 23:10:41 -0800
Fix Hysteresis Simulation
Diffstat:
5 files changed, 98 insertions(+), 95 deletions(-)
diff --git a/Paper/420_paper.pdf b/Paper/420_paper.pdf
Binary files differ.
diff --git a/Paper/420_paper.tex b/Paper/420_paper.tex
@@ -119,7 +119,7 @@ To quote Brian Eno, "Whatever you now find weird, ugly, uncomfortable, and nasty
about a new medium will surely become its signature. CD distortion, the jitteriness
of digital video, the crap sound of 8-bit - all of these will be cherished
and emulated as soon as they can be avoided." \cite{Eno}. The following describes
-a physical model of an analog tape machine, starting from basic physical
+a physical model of a Sony TC-260 tape machine, starting from basic physical
principles.
\section{Continuous Time System}
@@ -262,10 +262,10 @@ or,
\subsection{Hysteresis}
Beginning from \cref{eq5}, we can find the derivative of $M$ w.r.t. time,
as in \cite{Hysteresis}:
-\begin{multline}
- \frac{dM}{dt} = \frac{(1-c) \delta_M (M_sL(Q) - M)}{(1-c) \delta k - \alpha (M_sL(Q) - M)} \frac{dH}{dt} + \\
- c \frac{M_s}{a} \Big(\frac{dH}{dt} + \alpha \frac{dM}{dt} \Big) L'(Q)
-\end{multline}
+\begin{equation}
+ \frac{dM}{dt} = \frac{\frac{(1-c) \delta_M (M_sL(Q) - M)}{(1-c) \delta k - \alpha (M_sL(Q) - M)} \frac{dH}{dt} + c \frac{M_s}{a} \frac{dH}{dt} L'(Q)}{1 - c \alpha \frac{M_s}{a} L'(Q)}
+ \label{eq:dmdt}
+\end{equation}
%
where $Q = \frac{H + \alpha M}{a}$, and
@@ -273,52 +273,65 @@ where $Q = \frac{H + \alpha M}{a}$, and
L'(x) = \frac{1}{x^2} - \coth^2(x) + 1
\end{equation}
%
-Taking the second derivative, we find:
-\begin{multline}
- \frac{d^2 M}{dt^2} = \frac{(1-c) \delta_M (M_sL(Q) - M)}{(1-c) \delta k - \alpha (M_sL(Q) - M)} \frac{d^2H}{dt^2}\\
- + \frac{(1-c) \delta_M (c \frac{M_s}{a} (\frac{dH}{dt} + \alpha \frac{dM}{dt}) L'(Q) - \dot{M})}{(1-c) \delta k - \alpha (M_sL(Q) - M)} \frac{dH}{dt}\\
- + \frac{(1-c) \delta_M (M_sL(Q) - M)(-\alpha (c \frac{M_s}{a} (\frac{dH}{dt} + \alpha \frac{dM}{dt}) L'(Q) - \dot{M}))}{((1-c) \delta k - \alpha (M_sL(Q) - M))^2} \frac{dH}{dt}\\
- + c \frac{M_s}{a} \Big(\frac{dH}{dt} c \frac{M_s}{a} (\frac{dH}{dt} + \alpha \frac{dM}{dt}) L''(Q) + \frac{d^2H}{dt^2} L'(Q) + \alpha \frac{dM}{dt} c \frac{M_s}{a} (\frac{dH}{dt} + \alpha \frac{dM}{dt}) L''(Q) + \alpha \frac{d^2M}{dt^2} L'(Q) \Big)
-\end{multline}
-%
-where
-\begin{equation}
- L''(x) = 2 \coth(x) \cdot (\coth^2(x) - 1) - \frac{2}{x^3}
-\end{equation}
-%
-Solving for $\frac{dM}{dt}$ and $\frac{d^2M}{dt^2}$, we get:
+Note that \cref{eq:dmdt} can also be written in the general form for non-linear
+Ordinary Differential Equations:
\begin{equation}
- \frac{dM}{dt} = \frac{\frac{(1-c) \delta_M (M_sL(Q) - M)}{(1-c) \delta k - \alpha (M_sL(Q) - M)} \frac{dH}{dt} + c \frac{M_s}{a} \frac{dH}{dt} L'(Q)}{1 - c \alpha \frac{M_s}{a} L'(Q)} = f(t, M, \vec{u})
+ \frac{dM}{dt} = f(t,M,\vec{u})
\end{equation}
-\begin{multline}
- \frac{d^2 M}{dt^2} = \frac{1}{1 - c \alpha \frac{M_s}{a} L'(Q)} \Bigg[
-\frac{(1-c) \delta_M (M_sL(Q) - M)}{(1-c) \delta k - \alpha (M_sL(Q) - M)} \frac{d^2H}{dt^2}\\
-+ \frac{(1-c) \delta_M (c \frac{M_s}{a} (\frac{dH}{dt} + \alpha \frac{dM}{dt}) L'(Q) - \dot{M})}{(1-c) \delta k - \alpha (M_sL(Q) - M)} \frac{dH}{dt}\\
-+ \frac{(1-c) \delta_M (M_sL(Q) - M)(-\alpha (c \frac{M_s}{a} (\frac{dH}{dt} + \alpha \frac{dM}{dt}) L'(Q) - \dot{M}))}{((1-c) \delta k - \alpha (M_sL(Q) - M))^2} \frac{dH}{dt}\\
-+ c \frac{M_s}{a} \Big(\frac{d^2H}{dt^2} L'(Q) + c \frac{M_s}{a} L''(Q) (\dot{H} + \alpha \dot{M})^2 \Big) \Bigg]
-\end{multline}
-%
where $\vec{u} = \begin{bmatrix}
H \\
\dot{H}
\end{bmatrix}$.
\newline\newline
-Using the trapezoidal rule:
+Using the trapezoidal rule for derivative approximation, we find:
\begin{equation}
- \dot{\hat{H}}(n) = 2 \frac{\hat{H}(n) - \hat{H}(n-1)}{T} - \dot{\hat{H}}(n-1)
+ \dot{\hat{H}}(n) = 2\frac{\hat{H}(n) - \hat{H}(n-1)}{T} - \dot{\hat{H}}(n-1)
\end{equation}
+%
+We can use the Runge-Kutta 4th order method \cite{Yeh} to find an explicit solution
+for $\hat{M}(n)$:
+\begin{align}
+\begin{split}
+ k_1 &= T f \Big(n-1, \hat{M}(n-1), \hat{\vec{u}}(n-1) \Big)\\
+ k_2 &= T f \Big(n - \frac{1}{2}, \hat{M}(n-1) + \frac{k_1}{2}, \hat{\vec{u}} \Big(n-\frac{1}{2} \Big) \Big)\\
+ k_3 &= T f \Big(n- \frac{1}{2}, \hat{M}(n-1) + \frac{k_2}{2}, \hat{\vec{u}} \Big(n-\frac{1}{2} \Big) \Big)\\
+ k_4 &= T f \Big(n, \hat{M}(n-1) + k_3, \hat{\vec{u}}(n) \Big)\\
+ \hat{M}(n) &= \hat{M}(n-1) + \frac{k_1}{6} + \frac{k_2}{3} + \frac{k_3}{3} + \frac{k_4}{6}
+\end{split}
+\end{align}
+%
+We use linear interpolation to find the half-sample values used to calculate $k_2$ and $k_3$.
+
+\subsubsection{Numerical Considerations}
+To account for rounding errors in the Langevin function for values close to
+zero, we use the following for numerical evaluation, as in \cite{Hysteresis}:
\begin{equation}
- \ddot{\hat{H}}(n) = 2 \frac{\dot{\hat{H}}(n) - \dot{\hat{H}}(n-1)}{T} - \ddot{\hat{H}}(n-1)
+ L(x) = \begin{cases}
+ \coth(x) - \frac{1}{x} & \text{for $|x| > 10^{-4}$} \\
+ \frac{x}{3} & \text{otherwise}
+ \end{cases}
+\end{equation}
+\begin{equation}
+ L'(x) = \begin{cases}
+ \frac{1}{x^2} - \coth^{2}(x) + 1 & \text{for $|x| > 10^{-4}$} \\
+ \frac{1}{3} & \text{otherwise}
+ \end{cases}
\end{equation}
-%
-Then, using the semi-implicit trapezoidal rule for non-linear differential
-equations \cite{Yeh}.
-\begin{multline}
- \hat{M}(n) = \hat{M}(n-1)\\
- + \frac{T}{2} \frac{f[n, \hat{M}(n-1), \hat{\vec{u}}(n)] + f[n-1, \hat{M}(n-1), \hat{\vec{u}}(n-1)]}{1 - \frac{T}{2}\ddot{\hat{M}}(n-1)}
-\end{multline}
+\subsubsection{Simulation}
+The digitized hysteresis loop was implemented and tested offline
+in \texttt{Python}, using the constants $M_s$, $a$, $\alpha$, $k$,
+and $c$ from \cite{JilesAtherton1986}. For a sinusoidal input signal
+with frequency 2kHz, and varying amplitude from 800 - 2000 Amperes per
+meter, the following plot shows the Magnetisation output.
+
+\begin{figure}[ht]
+ \center
+ \includegraphics[width=3in]{../Simulations/Sim2-M_H.png}
+ \caption{\label{flowchart}{\it Digitized Hysteresis Loop}}
+\end{figure}
+%
\subsection{Play Head}
By combining \cref{eq11} with \cref{eq12,eq15}, we get:
\begin{equation}
diff --git a/Paper/references.bib b/Paper/references.bib
@@ -78,3 +78,16 @@ booktitle = {Theory of Magnetic Recording, by H.~Neal Bertram, pp.~372.~ISBN 052
year = 2009,
month = 6
}
+
+@ARTICLE{JilesAtherton1986,
+ author = {{Jiles}, D.~C. and {Atherton}, D.~L.},
+ title = "{Theory of ferromagnetic hysteresis}",
+ journal = {Journal of Magnetism and Magnetic Materials},
+ year = 1986,
+ month = sep,
+ volume = 61,
+ pages = {48-60},
+ doi = {10.1016/0304-8853(86)90066-1},
+ adsurl = {http://adsabs.harvard.edu/abs/1986JMMM...61...48J},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
diff --git a/Simulations/Sim2-M_H.png b/Simulations/Sim2-M_H.png
Binary files differ.
diff --git a/Simulations/hystersis.py b/Simulations/hystersis.py
@@ -1,12 +1,15 @@
import numpy as np
import matplotlib.pyplot as plt
+from matplotlib import animation
-T = 48000
-M_s = 1.6e6
-a = 1.1e3
+fs = 48000
+T = 1/fs # sample interval
+M_s = 350000 # Jiles book
+a = 1.1e3 # adjustable parameter
alpha = 1.6e-3
k = 4.0e2
c = 1.7e-1
+#H_c (coercivity) = 25000 - 30000 A/m (Jiles book)
# Langevin function
def L (x):
@@ -22,16 +25,9 @@ def L_d (x):
else:
return (1 / 3)
-# Langevin 2nd derivative
-def L_d2 (x):
- if (abs(x) > 10 ** -3):
- return 2 * (1 / np.tanh (x)) * ((1 / np.tanh (x)) ** 2 - 1) - (2 / x ** 3)
- else:
- return -(2 / 15) * x
-
# trapezoidal rule derivative
def deriv (x_n, x_n1, xDeriv_n1):
- return (2 / T) * (x_n - x_n1) - xDeriv_n1
+ return ((2 / T) * (x_n - x_n1)) - xDeriv_n1
# dM/dt or "non-linear function"
def f (M, H, H_d):
@@ -51,45 +47,22 @@ def f (M, H, H_d):
return (t1 + t2) / denominator
-# d^2M/dt^2
-def M_d2 (M, M_d, H, H_d, H_d2):
- Q = (H + alpha * M) / a
- M_diff = M_s * L (Q) - M
- delta = 1 if H_d > 0 else -1
- delta_M = 1 if np.sign (delta) == np.sign (M_diff) else 0
- L_prime = L_d (Q)
-
- denominator = 1 - c * alpha * (M_s / a) * L_prime
-
- t1_num = (1 - c) * delta_M * M_diff
- t1_den = (1 - c) * delta * k - alpha * M_diff
- t1 = (t1_num / t1_den) * H_d2
-
- M_diff2 = c * (M_s / a) * (H_d + alpha * M_d) * L_d (Q) - M_d
- t2_num = (1 - c) * delta_M * M_diff2
- t2 = (t2_num / t1_den) * H_d
-
- t3_num = (1 - c) * delta_M * M_diff * (-alpha) * M_diff2
- t3 = (t3_num / t1_den**2) * H_d
-
- L_2prime = L_d2 (Q)
- t4 = c * (M_s / a) * (H_d2 * L_prime + c * (M_s / a) * L_2prime * (H_d + alpha * M_d)**2)
-
- return (t1 + t2 + t3 + t4) / denominator
-
-def M_n (M_n1, f_n, f_n1, M_d2_n1):
- return M_n1 + (T / 2) * (f_n + f_n1) / (1 - (T/2) * M_d2_n1)
-
-t = np.linspace (0, 10, T * 50)
-H_in = np.concatenate ((1000 * np.sin (2 * np.pi * 2 * t[0:T*5]), 900 * np.sin (2 * np.pi * 2 * t[T*5:T*10]), \
- 800 * np.sin (2 * np.pi * 2 * t[T*10:T*15]), 750 * np.sin (2 * np.pi * 2 * t[T*15:T*20]), \
- 700 * np.sin (2 * np.pi * 2 * t[T*20:T*25]), 650 * np.sin (2 * np.pi * 2 * t[T*25:T*30]), \
- 600 * np.sin (2 * np.pi * 2 * t[T*30:T*35]), 550 * np.sin (2 * np.pi * 2 * t[T*35:T*40]),
- 500 * np.sin (2 * np.pi * 2 * t[T*40:T*45]), 400 * np.sin (2 * np.pi * 2 * t[T*45:T*50])))
+def M_n (M_n1, k1, k2, k3, k4):
+ return M_n1 + (k1 / 6) + (k2 / 3) + (k3 / 3) + (k4 / 6)
+
+#input signal
+t = np.linspace (0, 1, fs * 50)
+#H_in = 100 * np.sin (2 * np.pi * 100 * t)
+freq = 2000
+H_in = np.concatenate ((2000 * np.sin (2 * np.pi * freq * t[0:fs*5]), 1800 * np.sin (2 * np.pi * freq * t[fs*5:fs*10]), \
+ 1600 * np.sin (2 * np.pi * freq * t[fs*10:fs*15]), 1500 * np.sin (2 * np.pi * freq * t[fs*15:fs*20]), \
+ 1400 * np.sin (2 * np.pi * freq * t[fs*20:fs*25]), 1300 * np.sin (2 * np.pi * freq * t[fs*25:fs*30]), \
+ 1200 * np.sin (2 * np.pi * freq * t[fs*30:fs*35]), 1100 * np.sin (2 * np.pi * freq * t[fs*35:fs*40]), \
+ 1000 * np.sin (2 * np.pi * freq * t[fs*40:fs*45]), 800 * np.sin (2 * np.pi * freq * t[fs*45:fs*50])))
# plt.plot (t, H_in)
# plt.show()
-M_out = np.zeros (T * 50)
+M_out = np.zeros (fs * 50)
M_n1 = 0
H_n1 = 0
H_d_n1 = 0
@@ -101,29 +74,33 @@ for H in H_in:
H_d = deriv (H, H_n1, H_d_n1)
H_d2 = deriv (H_d, H_d_n1, H_d2_n1)
- f_n = f (M_n1, H, H_d)
- f_n1 = f (M_n1, H_n1, H_d_n1)
- M_d2_n1 = M_d2 (M_n1, f_n1, H_n1, H_d_n1, H_d2_n1)
+ k1 = T * f (M_n1, H_n1, H_d_n1)
+ k2 = T * f (M_n1 + k1/2, (H + H_n1) / 2, (H_d + H_d_n1) / 2)
+ k3 = T * f (M_n1 + k2/2, (H + H_n1) / 2, (H_d + H_d_n1) / 2)
+ k4 = T * f (M_n1 + k3, H, H_d)
- M = M_n (M_n1, f_n, f_n1, M_d2_n1)
+ M = M_n (M_n1, k1, k2, k3, k4)
M_n1 = M
H_n1 = H
H_d_n1 = H_d
H_d2_n1 - H_d2
- M_out[n] = (M)
+ M_out[n] = M
n += 1
- curPercent = (int) (n / (T * 50) * 100)
- if (curPercent > percent):
+ curPercent = (int) (n / (fs * 50) * 100)
+ if (curPercent > percent + 4):
percent = curPercent
print (str (percent) + "% completed")
MH = plt.figure (1)
-plt.plot (H_in, M_out)
-#MH.show()
+plt.plot (H_in / 1000, M_out / M_s)
+plt.xlabel ("Magnetic Field (A/m)")
+plt.ylabel ("Tape Magnetisation (A/m)")
+plt.title ("Simulated Ditigal Tape Magnetization Hysteresis Loop")
+MH.show()
Mt = plt.figure (2)
-plt.plot (t, M_out)
+plt.plot (t, M_out / M_s)
plt.show()
