# Skin Effect Load Update

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:

where

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):

with

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]):

which matches the definition of \(q\) in [2] (see above).

for \(Z_\Int\) of the Fortran implementation we get:

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:

From Wikipedia:

This is a factor of

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):

With the special cases of order 0 and 1:

The Wikipedia formula passes \(kr\) to the Bessel functions:

while the NEC-2 passes

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.