---
jupyter:
  description: 'The skin effect with planex2d and NGSolve: alternating current in a round copper wire against the
    Bessel solution and its resistance ratio from 50 Hz to 100 kHz, and a flat bar of the same area, where the current
    gathers at the edges.'
---

# The skin effect

Why does an alternating current flow near the surface of a wire? The current's changing magnetic field induces eddy
currents in the conductor, which oppose the current in the middle and add to it at the surface. The current density
decays inward over the skin depth $\delta = \sqrt{2/(\omega\mu\sigma)}$, which for copper is 9.3 mm at 50 Hz and
0.21 mm at 100 kHz. In a round wire of radius $a$, the current density follows a Bessel function, $J(r) \propto
J_0(kr)$ with $k = (1 - i)/\delta$, and the resistance rises over the direct current's by

$$\frac{R_{ac}}{R_{dc}} = \operatorname{Re} \frac{ka\,J_0(ka)}{2\,J_1(ka)}.$$

Here, NGSolve solves for the current in a copper wire of 1 mm radius and checks both against the Bessel solution. A
flat bar of nearly the same area follows, for which no closed form exists.

![The current density in a copper wire of 1 mm radius at 1, 10, 30 and 100 kHz, one quarter each](../_static/gallery/skin-effect.png){.result}

```python
import math

import matplotlib.pyplot as plt
import ngsolve as ngs
import numpy as np
import planex as px
from scipy.special import jv
```

```python nbsphinx="hidden"
from docs_figures import MAP, Chart, Row


def sampled(J, mesh, label):
    """|J| in A/mm² for Row, at points inside the conductor"""
    def field(face, plot=True, **_):
        (u0, v0), (u1, v1) = face.bounding_box
        gu, gv = np.meshgrid(np.linspace(u0, u1, 351), np.linspace(v0, v1, 91))
        uv = np.c_[gu.ravel(), gv.ravel()]
        uv = uv[face.Distance(uv, signed=True) < -1e-6]
        return {"points": [tuple(p) for p in uv], "values": np.abs(J(mesh(uv[:, 0], uv[:, 1]))).ravel() / 1e6,
                "label": label}
    return field
```

## The equations

In the cross-section, the magnetic vector potential $A$ along the conductor satisfies
$-\nabla \cdot (\mu_0^{-1} \nabla A) + i\omega\sigma A = \sigma E$ in the copper, and the same without the right-hand
side and the $\sigma$ term in the air. $E = 1\,\text{V/m}$ is the field along the conductor that drives the current,
and the current density is $J = \sigma (E - i\omega A)$. The resistance per meter follows from the power that the
current dissipates, $R = \int |J|^2/\sigma\,dA \,/\, |I|^2$, and the direct current's is $1/(\sigma A_{cu})$ for the
copper's area $A_{cu}$.

`section` glues the conductor into a circle of air of radius 20 mm, with a fine mesh along the conductor's edges, and
`current` solves for one frequency:

```python
sigma, mu0 = 5.8e7, 4e-7 * math.pi


def section(conductor):
    conductor.name = "copper"
    conductor.edges.maxh = 2.5e-5
    air = px.Circle((0, 0), 0.02).Face()
    air.name = "air"
    air.edges.name = "outer"
    mesh = ngs.Mesh(px.OCCGeometry(px.Glue([air, conductor])).GenerateMesh(maxh=5e-4))
    mesh.Curve(3)
    return mesh


def current(mesh, f):
    omega = 2 * math.pi * f
    copper = mesh.Materials("copper")
    V = ngs.H1(mesh, order=3, complex=True, dirichlet="outer")
    u, v = V.TnT()
    K = ngs.BilinearForm(ngs.grad(u) * ngs.grad(v) / mu0 * ngs.dx
                         + 1j * omega * sigma * u * v * ngs.dx(definedon=copper)).Assemble()
    F = ngs.LinearForm(sigma * v * ngs.dx(definedon=copper)).Assemble()
    A = ngs.GridFunction(V)
    A.vec.data = K.mat.Inverse(V.FreeDofs(), inverse="umfpack") * F.vec
    J = sigma * (1 - 1j * omega * A)
    I = ngs.Integrate(J, mesh, definedon=copper)
    loss = ngs.Integrate(ngs.Norm(J) ** 2 / sigma, mesh, definedon=copper)
    area = ngs.Integrate(1, mesh, definedon=copper)
    return J, loss / abs(I) ** 2 * sigma * area
```

## A round wire

The wire has a radius of 1 mm. From 50 Hz to 100 kHz, the skin depth falls from more than nine times the radius to a
fifth of it, and the resistance ratio follows the Bessel formula:

```python
r0 = 1e-3
wire = section(px.Circle((0, 0), r0).Face())
profiles = {}
for f in (50, 1e3, 1e4, 1e5):
    J, ratio = current(wire, f)
    k = (1 - 1j) * math.sqrt(math.pi * f * mu0 * sigma)
    exact = (k * r0 * jv(0, k * r0) / (2 * jv(1, k * r0))).real
    profiles[f] = (J, k)
    print("f", f, "Hz | skin depth", round(math.sqrt(2) / abs(k) * 1e3, 3), "mm | R_ac/R_dc", ratio, "exact", exact)
```

Across the wire, the current density at 50 Hz is nearly uniform. At 100 kHz it falls steeply from the surface, to about
an eighth of its value halfway to the center, and the computed densities lie on the Bessel solution (Figure 1):

```python nbsphinx="hide-input"
r, few = np.linspace(0, 0.999 * r0, 200), np.linspace(0, 0.999 * r0, 12)
surface = {f: abs(J(wire(0.999 * r0, 0.0))) for f, (J, k) in profiles.items()}
Chart("radius in mm", "|J| against the surface",
      dots={f"finite elements, {f:g} Hz": (few * 1e3, np.abs(J(wire(few, 0 * few))).ravel() / surface[f])
            for f, (J, k) in profiles.items()},
      lines={f"Bessel, {f:g} Hz": (r * 1e3, np.abs(jv(0, k * r) / jv(0, k * 0.999 * r0)))
             for f, (J, k) in profiles.items()},
      title="The current density across a copper wire of 1 mm radius")
```

## A flat bar

A flat bar of 3.5 × 0.9 mm has nearly the wire's area, 3.15 mm² against 3.14 mm², but almost half again its
perimeter, 8.8 mm against 6.3 mm. At 100 kHz, its current gathers at the two narrow ends and most of all at the
corners (Figure 2), and its resistance rises less:

```python
flat = px.QuickSketch().RectangleC(3.5e-3, 0.9e-3).Face()
bar = section(flat)
J, ratio = current(bar, 1e5)
print("flat bar at 100 kHz: R_ac/R_dc", ratio, "| round wire", current(wire, 1e5)[1])
```

```python nbsphinx="hide-input"
Row((flat, "The current density at 100 kHz", {"field": sampled(J, bar, "|J| in A/mm²"), "cmap": MAP.reversed(),
                                               "bar": "right"}))
```

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
from gallery import save_figure
from docs_figures import INK


quarters = [current(wire, f)[0] for f in (1e3, 1e4, 3e4, 1e5)]


def quadrants(face, plot=True, **_):
    """|J| against the surface for PlotField, one frequency per quarter of the wire, rising counterclockwise"""
    g = np.linspace(-r0, r0, 301)
    uv = np.array([(x, y) for y in g for x in g])
    uv = uv[np.hypot(uv[:, 0], uv[:, 1]) < 0.999 * r0]
    quarter = (np.arctan2(uv[:, 1], uv[:, 0]) % (2 * np.pi) // (np.pi / 2)).astype(int)
    values = np.empty(len(uv))
    for q, J in enumerate(quarters):
        at = uv[quarter == q]
        values[quarter == q] = np.abs(J(wire(at[:, 0], at[:, 1]))).ravel() / abs(J(wire(0.999 * r0, 0.0)))
    return {"points": [tuple(p) for p in uv], "values": values, "label": "|J| against the surface"}


fig, ax = plt.subplots(figsize=(5, 4))
copper = px.Circle((0, 0), r0).Face()
px.PlotField(quadrants, copper, ax=ax, colormap=MAP.reversed(), min=0, max=1)
px.Plot(px.Circle((0, 0), 1.35 * r0).Face() - copper, ax=ax, face_col="#c3cbd3", line_col=INK, arc_col=INK)
ax.plot([-r0, r0], [0, 0], color="white", lw=1.5)
ax.plot([0, 0], [-r0, r0], color="white", lw=1.5)
ax.set_aspect("equal", adjustable="box")
ax.set_xlim(-1.4 * r0, 1.4 * r0)
ax.set_ylim(-1.4 * r0, 1.4 * r0)
save_figure(fig, "skin-effect", bare=True, bars={1: "|J| against the surface"})
```