.. title: Skin Effect Load Update
.. slug: 28
.. date: 2024-07-28 20:45
.. tags: documentation,english,howto,hamradio,antenna modeling
.. description:
.. wp-status: publish
.. has_math: true
.. |--| unicode:: U+2013 .. en dash
.. |ohm| unicode:: U+02126 .. Omega
:trim:
.. |_| unicode:: U+00A0 .. Non-breaking space
:trim:
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:
.. math::
\newcommand{\Int}{{\mathrm\scriptscriptstyle int}}
\newcommand{\ber}{\mathop{\mathrm{ber}}\nolimits}
\newcommand{\bei}{\mathop{\mathrm{bei}}\nolimits}
.. math::
Z_\Int = \frac{k\rho}{2\pi r}\frac{J_0 (kr)}{J_1 (kr)}
where
.. math::
k = \sqrt{\frac{-j\omega\mu}{\rho}}
and :math:`r` is the wire radius, :math:`J_v` are the Bessel functions of
the first kind of order :math:`v`. :math:`Z_\Int` is the impedance *per
unit length* of wire.
`NEC-2`_ uses a different formulation [2]_ (p.75):
.. math::
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]
with
.. math::
q = (\omega\mu\sigma)^{1/2} a_i
and :math:`a_i` the wire radius, :math:`\Delta_i` the length of the wire
segment, :math:`\sigma` the wire conductivity and :math:`\ber` and
:math:`\bei` the Kelvin functions. In the following we will use
:math:`r` for the wire radius, not :math:`a_i`. Note that
:math:`\ber^\prime` is the derivative of :math:`\ber`, we can also write
:math:`\ber` as :math:`\ber_0` (order 0) and :math:`\ber^\prime` as
:math:`\ber_1` (order 1) and similar for :math:`\bei`.
In the following I'll be a little bit sloppy about the usage of the
permeability of free space :math:`\mu_0` vs. the permeability of the
wire :math:`\mu`, for non-ferromagnetic wires these are really (almost)
the same: :math:`\mu = \mu_r * \mu_0` with :math:`\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
:math:`\lambda` (parameter ``SIGL``) and the radius divided by the
wavelength :math:`\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 :math:`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`` :math:`\approx \frac{\mu_0 c}{2\pi}`
``TPCMU`` :math:`\approx 2\pi c\mu_0`
(:math:`c` is the speed of light) we get for ``X`` from the code (which
we call :math:`q` now consistent with the formula from [2]_):
.. math::
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} \\
which matches the definition of :math:`q` in [2]_ (see above).
for :math:`Z_\Int` of the Fortran implementation we get:
.. math::
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] \\
Which does *NOT* match the formula for :math:`Z_i` from [2]_ above.
It looks like the formula from [2]_ (p. 75) should have the frequency
:math:`f` instead of :math:`\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 :math:`\Delta_i`, making it
a *per unit length* formula like the one from Wikipedia), also note that
:math:`\rho=1/\sigma`:
From code:
.. math::
\frac{j}{2\pi r}\sqrt{\frac{\omega\mu}{\sigma}}
From Wikipedia:
.. math::
\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}}
This is a factor of
.. math::
\frac{\sqrt{-j}}{j} = -\frac{1+j}{\sqrt{2}}
because :math:`\sqrt{-j} = \frac{1-j}{\sqrt{2}}`
For the Kelvin functions we have the following equivalence to Bessel
functions [3]_ (10.61, `Kelvin Functions`_):
.. math::
\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)
With the special cases of order 0 and 1:
.. math::
\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) \\
The Wikipedia formula passes :math:`kr` to the Bessel functions:
.. math::
kr = r\sqrt{\frac{-j\omega\mu}{\rho}} = r\sqrt{-j\omega\mu\sigma}
while the `NEC-2`_ passes
.. math::
q = r\sqrt{\omega\mu\sigma}
which differ by the factor :math:`\sqrt{-j}`. This means that we pass a
factor of :math:`\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 :math:`\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.
.. iframe:: /content/2024/07/plotly/z_resist_2.html
.. iframe:: /content/2024/07/plotly/z_react_2.html
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.
.. iframe:: /content/2024/07/plotly/resistivity2.html
.. [1] Ralf Schlatterbeck. `Skin effect load and impedance`_. Blog post,
Open Source Consulting, July 2024.
.. [2] G. J. Burke and A. J. Poggio. Numerical electromagnetics code (NEC)
|--| method of moments, `Part I: Program description`_ |--| theory.
January 1981.
.. [3] `Digital library of mathematical functions`_. NIST, June 2024.
.. _pymininec: https://github.com/schlatterbeck/pymininec
.. _Mininec: https://github.com/schlatterbeck/MiniNec
.. _`NEC-2`: https://en.wikipedia.org/wiki/Numerical_Electromagnetics_Code
.. _`skin effect`: https://en.wikipedia.org/wiki/Skin_effect
.. _`Skin effect load and impedance`: https://blog.runtux.com/posts/2024/07/04/
.. _`Part I: Program description`: https://nec2.org/other/nec2prt1.pdf
.. _`Digital library of mathematical functions`: https://dlmf.nist.gov/
.. _`Kelvin Functions`: https://dlmf.nist.gov/10.61
.. _WINE: https://www.winehq.org/