Sharp or blunt lightning rods

Why are lightning rods pointed? At the surface of a conductor, the electric field is strongest where the surface is curved most sharply, and Benjamin Franklin’s rods ended in sharp points. For a rod shaped as a half ellipse of height \(a\) and half width \(b\) on flat ground, the field at the tip has a closed form. With its mirror image below the ground, the rod is a conducting elliptic cylinder in a uniform field \(E_0\), and in elliptic coordinates \(\xi, \eta\) its potential is \(\varphi = E_0 c \cos\eta \,(\cosh\xi - \cosh\xi_0\, e^{\xi_0 - \xi})\), with \(c = \sqrt{a^2 - b^2}\) and the rod’s surface at \(\xi_0\). At the tip, this gives

\[E_{tip} = E_0 \left(1 + \frac{a}{b}\right),\]

so the sharper the rod, the stronger the field at its tip, without a limit. Field tests tell a different story. In seven years on a mountain, Moore, Rison, Mathis and Aulich (2000) saw twelve blunt rods struck and none of the sharp rods beside them, and their calculations showed why: a few centimeters above the tip, a blunt rod’s field is stronger than a sharper rod’s of the same height.

Here, NGSolve solves the electric potential around four rods of the same height, from a half circle to a sharp half ellipse. The tip fields are checked against \(1 + a/b\), the fields above the tips are compared, and a rod goes onto the ridge of a house.

The electric field around a blunt and a sharp lightning rod of the same height

[1]:
import math

import matplotlib.pyplot as plt
import ngsolve as ngs
import numpy as np
import planex as px

The rods

Every rod is 10 high, and the half widths run from 10, a half circle, down to 1.25. px.Ellipse draws each rod exactly, standing upright with angle=90, and the ground cuts it in half. The air is a box 600 high and 1200 wide, large against the rods. The rod’s edges and the box’s ground and top are named before the Boolean operation, which keeps the names. The rod’s edges get a mesh size of a quarter of the radius of curvature at the tip, \(b^2/a\):

[3]:
a, L = 10.0, 600.0


def air(b):
    rod = px.Face([px.Ellipse((0, 0), a, b, angle=90)])
    rod.edges.name = "rod"
    rod.edges.maxh = b**2 / a / 4
    box = px.QuickSketch((-L, 0)).Rectangle(2 * L, L).Face()
    box.edges(near=(0, 0)).name = "ground"
    box.edges(near=(0, L)).name = "top"
    return (box - rod).Face()

The field

The rod and the ground are at the potential 0, and the top of the box at \(E_0 L\) with \(E_0 = 1\), so the field far from the rod is uniform. NGSolve solves Laplace’s equation for the potential with elements of order 3 on a mesh curved to the ellipse, and field returns the strength of the field at given points:

[4]:
def field(b, around=air):
    mesh = ngs.Mesh(px.OCCGeometry(around(b)).GenerateMesh(maxh=L / 10, grading=0.25))
    mesh.Curve(3)
    V = ngs.H1(mesh, order=3, dirichlet="rod|ground|top")
    u, v = V.TnT()
    phi = ngs.GridFunction(V)
    phi.Set(L, definedon=mesh.Boundaries("top"))
    K = ngs.BilinearForm(ngs.grad(u) * ngs.grad(v) * ngs.dx).Assemble()
    phi.vec.data += K.mat.Inverse(V.FreeDofs(), inverse="umfpack") * (-K.mat * phi.vec)
    strength = ngs.Norm(ngs.grad(phi))
    return lambda h, v: strength(mesh(h, v))


fields = {b: field(b) for b in (10, 5, 2.5, 1.25)}
for b, E in fields.items():
    print("b", b, "| tip field", E(0.0, a + 1e-6), "| exact", 1 + a / b)
b 10 | tip field 2.0014305040388076 | exact 2.0
b 5 | tip field 3.0005010923854116 | exact 3.0
b 2.5 | tip field 5.001500184658622 | exact 5.0
b 1.25 | tip field 9.000760132212367 | exact 9.0

Around the blunt and the sharp rod, the field is strongest at the tip (Figure 1). The color scale stops at 4, below the sharp rod’s 9:

Half width 10
(a) Half width 10
Half width 1.25
(b) Half width 1.25
Figure 1: The field around a blunt and a sharp rod of the same height

Above the tip

Along the axis above the rod, the exact solution gives the field as well, with \(y = c\cosh\xi\): \(E/E_0 = 1 + \cosh\xi_0\, e^{\xi_0 - \xi} / \sinh\xi\), and for the half circle \(1 + a^2/y^2\). The sharp rod’s field is the highest at its tip and falls fastest. A third of the rod’s height above the tip, the half circle’s field has overtaken it:

[6]:
def exact(y, b):
    if b == a:
        return 1 + a**2 / y**2
    c = math.sqrt(a**2 - b**2)
    xi0, xi = math.acosh(a / c), math.acosh(y / c)
    return 1 + math.cosh(xi0) * math.exp(xi0 - xi) / math.sinh(xi)


above = np.linspace(a + 0.01, 3 * a, 2000)
blunt, sharp = (fields[b](0 * above, above).ravel() for b in (10, 1.25))
crossing = above[np.argmax(blunt > sharp)] - a
exact_crossing = above[np.argmax([exact(y, 10) > exact(y, 1.25) for y in above])] - a
print("the half circle overtakes the sharp rod", round(crossing, 2), "above the tip, exact", round(exact_crossing, 2))
the half circle overtakes the sharp rod 3.41 above the tip, exact 3.41

The computed fields of the two rods follow the exact solution, and the crossing is marked (Figure 2):

field in units of E0 over height above the tip: finite elements, half width 10 0.01 2, 1.01 1.83, 2.01 1.69, 3.01 1.59, 4.01 1.51, 5.01 1.44, 6.01 1.39, 7.01 1.35, 8.01 1.31, 9.01 1.28, 10 1.25, 11 1.23, 12 1.21, 13 1.19, 14 1.17, 15 1.16, 16 1.15, 17 1.14, 18 1.13, 19 1.12; finite elements, half width 1.25 0.01 8.48, 1.01 2.49, 2.01 1.88, 3.01 1.62, 4.01 1.48, 5.01 1.38, 6.01 1.31, 7.01 1.26, 8.01 1.23, 9.01 1.2, 10 1.17, 11 1.15, 12 1.14, 13 1.12, 14 1.11, 15 1.1, 16 1.09, 17 1.09, 18 1.08, 19 1.07
Figure 2: The field on the axis above a blunt and a sharp rod

The sharp tip concentrates its field into a small region, and the blunt rod spreads its field farther into the air above it. This is the effect that Moore and his colleagues found in their calculations for round rods. The model here is planar, a ridge rather than a round rod, but the comparison comes out the same way.

On a roof

On a building, the rod stands on the ridge. A house 10 wide, with walls 6 high and a roof pitched at 35°, gets a rod 2 high and 0.2 in half width on its ridge. The house and the rod are one conductor at the potential 0, named rod before they are joined, and the same field solves the air around them:

[8]:
ridge = 6 + 5 * math.tan(math.radians(35))


def building(b, height=2.0):
    walls = px.QuickSketch((-5, 0)).Line((5, 0)).Line((5, 6)).Line((0, ridge)).Line((-5, 6)).Close().Face()
    walls.edges.name = "rod"
    rod = px.Face([px.Ellipse((0, ridge), height, b, angle=90)])
    rod.edges.name = "rod"
    rod.edges.maxh = b**2 / height / 4
    return (walls + rod).Face()


def house(b):
    box = px.QuickSketch((-L, 0)).Rectangle(2 * L, L).Face()
    box.edges(near=(0, 0)).name = "ground"
    box.edges(near=(0, L)).name = "top"
    return (box - building(b)).Face()


on_roof = field(0.2, house)
print("tip field on the ridge", on_roof(0.0, ridge + 2 + 1e-6), "| the same rod on flat ground", 1 + 2 / 0.2)
tip field on the ridge 22.555604205908654 | the same rod on flat ground 11.0

The house raises the field around its ridge, so the rod’s tip sees about twice the field of the same rod on flat ground (Figure 3). Away from the house, the field returns to \(E_0\):

A rod on the ridge
Figure 3: The field around a house with a lightning rod on its ridge