Skin Effect Load Update

Ralf Schlatterbeck

Also available in:

Deutsch

In my last post [1] 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 [2] (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 [2]):

\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 [2] (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 [2] above. It looks like the formula from [2] (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 [3] (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.