In meinem letzten Post habe ich mich gefragt, warum meine
Implementierung von Skin Effekt Belag (Drähte mit einem endlichen
Widerstand) nicht mit den Resultaten von NEC-2 übereinstimmten.
Es stellt sich heraus dass ich falsch lag in der Annahme dass
Skin-Effekt Belag nur eine rein ohmsche Komponente hat, durch den Skin
Effekt ergibt sich in Wirklichkeit eine relativ große reaktive
Komponente.
Hinweis: Dieser Artikel enthält ziemlich viel Mathematik, vielleicht
möchten Sie ja ans Ende scrollen und checken wie die neue
Implementierung in pymininec sich verglichen mit meinem letzten Artikel
verhält wo nur ein rein ohmscher Widerstand angenommen wurde.
Der englische Wikipedia Artikel zum skin effect enthält die
korrekte Formel:
\begin{equation*}
\newcommand{\Int}{{\mathrm\scriptscriptstyle int}}
\newcommand{\ber}{\mathop{\mathrm{ber}}\nolimits}
\newcommand{\bei}{\mathop{\mathrm{bei}}\nolimits}
\end{equation*}
\begin{equation*}
Z_\Int = \frac{k\rho}{2\pi r}\frac{J_0 (kr)}{J_1 (kr)}
\end{equation*}
wobei
\begin{equation*}
k = \sqrt{\frac{-j\omega\mu}{\rho}}
\end{equation*}
und \(r\) ist der Draht-Radius, \(J_v\) sind die
Bessel-Funktionen der ersten Gattung mit der Ordnung
\(v\). \(Z_\Int\) ist die Impedanz pro Längeneinheit des
Drahtes.
NEC-2 verwendet eine andere Formel (S.75):
\begin{equation*}
Z_i = \frac{j\Delta_i}{a_i} \sqrt{\frac{\omega\mu}{2\pi\sigma}}
\left[\frac{\ber (q) + j\bei (q)}
{\ber^\prime(q)+ j\bei^\prime(q)}\right]
\end{equation*}
mit
\begin{equation*}
q = (\omega\mu\sigma)^{1/2} a_i
\end{equation*}
und \(a_i\) ist der Draht-Radius, \(\Delta_i\) die Länge des
Drahtsegments, \(\sigma\) die Elektrische Leitfähigkeit des Drahtes
und \(\ber\) und \(\bei\) die Kelvin-Funktionen. Im folgenden
werde ich wieder \(r\) für den Draht-Radius, und nicht \(a_i\)
verwenden. Nachdem \(\ber^\prime\) die Ableitung von \(\ber\)
ist, können wir
\(\ber\) als \(\ber_0\) (Ordnung 0) und \(\ber^\prime\) als
\(\ber_1\) (Ordnung 1) schreiben und analog für \(\bei\).
Im folgenden bin ich etwas ungenau in der Verwendung der magnetischen
Permeabilität im Freiraum \(\mu_0\) vs. der Permeabilität eines
Drahtes \(\mu\), für nicht ferromagnetische Drähte sind diese (fast)
identisch: \(\mu = \mu_r * \mu_0\) mit \(\mu_r\approx 1\) für
nicht ferromagnetische Materialien.
Der NEC-2 code verwendet die Fortran Funktion ZINT mit zwei
Parametern, der Leitfähigkeit multipliziert mit der Wellenlänge
\(\lambda\) (Parameter SIGL) und der Radius geteilt durch die
Wellenlänge \(\lambda\) (Parameter ROLAM). Neben der
Approximation der Kelvin-Funktionen (der Term mit dem Bruch aus
Kelvin-Funktionen ist BR1), berechnet ZINT:
X = SQRT (TPCMU * SIGL) * ROLAM
ZINT = FJ * SQRT (CMOTP / SIGL) * BR1 / ROLAM
mit TPCMU = 2368.705, CMOTP=60.0, und FJ ist \(j\).
Das berechnete X wird der Berechnung der Kelvin-Funktionen übergeben
mit BR1 als Resultat des Bruch-Terms. Das berechnete ZINT
enthält nicht die Segment-Länge, ist also äquivalent mit der Wikipedia
Definition in pro Längeneinheit. Mit den Reverse-Engineerten magischen
Zahlen
CMOTP \(\approx \frac{\mu_0 c}{2\pi}\)
TPCMU \(\approx 2\pi c\mu_0\)
(\(c\) ist die Lichtgeschwindigkeit) bekommen wir für X aus dem
Code (was wir jetzt \(q\) nennen um mit der obigen Formel aus
konsistent zu sein):
\begin{align*}
q &= \frac{r}{\lambda}\sqrt{2368.705\sigma\lambda} \\
&= r\sqrt{\frac{2368.705\sigma}{\lambda}} \\
&\approx r\sqrt{\frac{2\pi c\mu_0\sigma}{\lambda}} \\
&= r\sqrt{\frac{2\pi c\mu_0\sigma f}{c}} \\
&= r\sqrt{\frac{2\pi \mu_0\sigma \omega}{2\pi}} \\
&= r\sqrt{\omega\mu_0\sigma} \\
\end{align*}
Was mit der obigen Definition von \(q\) in übereinstimmt.
Für \(Z_\Int\) der Fortran-Implementierung bekommen wir:
\begin{align*}
Z_\Int &= \frac{j}{r} \sqrt{\frac{\mu_0 c}{2\pi\sigma\lambda}}
\left[\frac{\ber_0 (q) + j\bei_0 (q)}
{\ber_1(q)+ j\bei_1(q)}\right] \\
&= \frac{j}{r} \sqrt{\frac{\mu_0 c f}{2\pi\sigma c}}
\left[\frac{\ber_0 (q) + j\bei_0 (q)}
{\ber_1(q)+ j\bei_1(q)}\right] \\
&= \frac{j}{r} \sqrt{\frac{\mu_0 f}{2\pi\sigma}}
\left[\frac{\ber_0 (q) + j\bei_0 (q)}
{\ber_1(q)+ j\bei_1(q)}\right] \\
&= \frac{j}{r} \sqrt{\frac{\mu_0 \omega}{2\pi 2\pi\sigma}}
\left[\frac{\ber_0 (q) + j\bei_0 (q)}
{\ber_1(q)+ j\bei_1(q)}\right] \\
&= \frac{j}{2\pi r} \sqrt{\frac{\mu_0 \omega}{\sigma}}
\left[\frac{\ber_0 (q) + j\bei_0 (q)}
{\ber_1(q)+ j\bei_1(q)}\right] \\
\end{align*}
Diese Formel stimmt NICHT mit der Formel für \(Z_i\) aus
weiter oben überein. Es sieht so aus als ob die Formel aus (S. 75)
die Frequenz \(f\) statt der Kreisfrequenz \(\omega\) unter der
Wurzel haben sollte. Die aus der NEC-2 Implementierung abgeleitete
Formel ist aber konsistent mit der obigen Formel aus Wikipedia.
Der Term vor dem Bruch mit den Kelvin-Funktionen unterscheidet sich um
einen Faktor von der Wikipedia Formel (Ich lasse \(\Delta_i\) weg
und habe so eine Formel pro Längeneinheit wie in Wikipedia). Ausserdem
gilt \(\rho=1/\sigma\):
Aus dem Code abgeleitet:
\begin{equation*}
\frac{j}{2\pi r}\sqrt{\frac{\omega\mu}{\sigma}}
\end{equation*}
Aus Wikipedia:
\begin{equation*}
\frac{k\rho}{2\pi r}
= \sqrt{\frac{-j\omega\mu}{\rho}}\frac{\rho}{2\pi r}
= \frac{1}{2\pi r}\sqrt{-j\omega\mu\rho}
= \frac{1}{2\pi r}\sqrt{\frac{-j\omega\mu}{\sigma}}
\end{equation*}
Das ist ein Faktor von
\begin{equation*}
\frac{\sqrt{-j}}{j} = -\frac{1+j}{\sqrt{2}}
\end{equation*}
weil \(\sqrt{-j} = \frac{1-j}{\sqrt{2}}\)
Für die Kelvin-Funktionen gilt die folgende Äquivalenz zu den Bessel
Funktionen (10.61, Kelvin Functions):
\begin{equation*}
\ber_v x + j \bei_v x = J_v \left(xe^{3\pi j/4}\right)
= e^{v\pi j} J_v \left(xe^{-\pi j / 4}\right)
\end{equation*}
Mit den Spezialfällen der Ordnung 0 und 1:
\begin{align*}
\ber_0 x + j \bei_0 x &= e^{0} J_0 \left(xe^{-\pi j / 4}\right) \\
&= J_0 \left(\frac{1-j}{\sqrt{2}}x\right) \\
\ber_1 x + j \bei_1 x &= e^{j\pi} J_1 \left(xe^{-\pi j / 4}\right) \\
&= -J_1 \left(\frac{1-j}{\sqrt{2}}x\right) \\
\end{align*}
Die Wikipedia Formel übergibt \(kr\) an die Bessel-Funktionen:
\begin{equation*}
kr = r\sqrt{\frac{-j\omega\mu}{\rho}} = r\sqrt{-j\omega\mu\sigma}
\end{equation*}
während NEC-2 folgendes übergibt:
\begin{equation*}
q = r\sqrt{\omega\mu\sigma}
\end{equation*}
Diese unterscheiden sich um den Faktor \(\sqrt{-j}\). Das bedeutet
wir übergeben einen Faktor um \(\sqrt{2}\) zu groß an die
Bessel-Funktionen. Daraus ergibt sich wegen der Symmetrie der Bessel
Funktionen und er Tatsache dass die Funktionen im Nenner die Ableitung
der Funktionen im Zähler sind, dass das Ergebnis um einen Faktor von
\(\sqrt{2}\) zu groß ist was den Faktor im ersten Term ausgleicht.
Nach der Implementierung der neuen Formel in pymininec verhält sich der
reaktive Teil der Impedanz ähnlich wie bei NEC-2 und der ohmsche
Widerstand ist noch näher an den Ergebnissen von NEC als vorher.
Auch wenn man den Imaginärteil (die Reaktanz) plottet sind die
Resultate jetzt wie erwartet. Die Werte von MMANA sind immer noch
überraschend, eine mögliche Ursache ist dass ich MMANA unter WINE in
Linux verwendet habe, die ohmschen Werte der Impedanz die von MMANA
berechnet wurden schauen allerdings korrekt aus.
