In my last post I was wondering why my implementation of skin effect
load in pymininec did not match the results we get from NEC-2. Turns
out I was wrong in assuming that the skin effect load is purely
resistive, it has a large reactive component.
Note: This article contains a lot of math. You may want to skip to the
end and check how the new implementation in pymininec behaves compared
to may last article which used only a purely resistive load from the
skin effect resistivity.
The english Wikipedia article on skin effect does have the correct
formula:
\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*}
where
\begin{equation*}
k = \sqrt{\frac{-j\omega\mu}{\rho}}
\end{equation*}
and \(r\) is the wire radius, \(J_v\) are the Bessel functions of
the first kind of order \(v\). \(Z_\Int\) is the impedance per
unit length of wire.
NEC-2 uses a different formulation (p.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*}
with
\begin{equation*}
q = (\omega\mu\sigma)^{1/2} a_i
\end{equation*}
and \(a_i\) the wire radius, \(\Delta_i\) the length of the wire
segment, \(\sigma\) the wire conductivity and \(\ber\) and
\(\bei\) the Kelvin functions. In the following we will use
\(r\) for the wire radius, not \(a_i\). Note that
\(\ber^\prime\) is the derivative of \(\ber\), we can also write
\(\ber\) as \(\ber_0\) (order 0) and \(\ber^\prime\) as
\(\ber_1\) (order 1) and similar for \(\bei\).
In the following I'll be a little bit sloppy about the usage of the
permeability of free space \(\mu_0\) vs. the permeability of the
wire \(\mu\), for non-ferromagnetic wires these are really (almost)
the same: \(\mu = \mu_r * \mu_0\) with \(\mu_r\approx 1\) for
non-ferromagnetic material.
The NEC-2 code, uses a Fortran function ZINT which gets
two parameters, the conductivity multiplied by the wavelength
\(\lambda\) (parameter SIGL) and the radius divided by the
wavelength \(\lambda\) (parameter ROLAM). Apart from an
approximation of the Kelvin functions (the term computing the fraction
of the Kelvin functions is BR1), ZINT computes:
X = SQRT (TPCMU * SIGL) * ROLAM
ZINT = FJ * SQRT (CMOTP / SIGL) * BR1 / ROLAM
with TPCMU = 2368.705, CMOTP=60.0, and FJ is \(j\).
The computed X is passed to the computation of the Kelvin functions
resulting in BR1 for the fraction term. Note that the computed
ZINT does not include the length, it is equivalent with the
Wikipedia definition in unit length.
With the reverse-engineered magic numbers
CMOTP \(\approx \frac{\mu_0 c}{2\pi}\)
TPCMU \(\approx 2\pi c\mu_0\)
(\(c\) is the speed of light) we get for X from the code (which
we call \(q\) now consistent with the formula from ):
\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*}
which matches the definition of \(q\) in (see above).
for \(Z_\Int\) of the Fortran implementation we get:
\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*}
Which does NOT match the formula for \(Z_i\) from above.
It looks like the formula from (p. 75) should have the frequency
\(f\) instead of \(\omega\) under the square root. The formula
derived from the NEC-2 implementation, however, is consistent with
the Wikipedia formula above.
The term before the fraction of Kelvin functions differs by a factor
from the Wikipedia formula (I'm leaving out \(\Delta_i\), making it
a per unit length formula like the one from Wikipedia), also note that
\(\rho=1/\sigma\):
From code:
\begin{equation*}
\frac{j}{2\pi r}\sqrt{\frac{\omega\mu}{\sigma}}
\end{equation*}
From 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*}
This is a factor of
\begin{equation*}
\frac{\sqrt{-j}}{j} = -\frac{1+j}{\sqrt{2}}
\end{equation*}
because \(\sqrt{-j} = \frac{1-j}{\sqrt{2}}\)
For the Kelvin functions we have the following equivalence to Bessel
functions (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*}
With the special cases of order 0 and 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*}
The Wikipedia formula passes \(kr\) to the Bessel functions:
\begin{equation*}
kr = r\sqrt{\frac{-j\omega\mu}{\rho}} = r\sqrt{-j\omega\mu\sigma}
\end{equation*}
while the NEC-2 passes
\begin{equation*}
q = r\sqrt{\omega\mu\sigma}
\end{equation*}
which differ by the factor \(\sqrt{-j}\). This means that we pass a
factor of \(\sqrt{2}\) too large to the Bessel functions. In turn
due to the symmetry and the fact that the functions in the denominator
are the derivative of the numerator, the output is too large by a factor
of \(\sqrt{2}\) which compensates the factor of the first term.
Implementing the new formula in pymininec now the reactance part of the
impedance behaves similar to NEC-2 and the results for the resistance are
even closer than before.
Also when plotting reactance by wire resistivity the results are now as
expected. The values for MMANA are still surprising, the reason might
be that I was running MMANA under WINE in Linux, the resistance values
computed by MMANA look OK, though.
