The Halbach array

Can a row of magnets have its field on one side only? In a Halbach array, the magnetization turns from block to block, by a quarter turn in the common version. On one side of the array the fields of the blocks add up, and on the other side they cancel. An ideal array, whose magnetization turns continuously along it, has no field at all on its weak side. On its strong side, the field has the constant magnitude \(B = B_r\,(1 - e^{-kd})\,e^{-ky}\), for the remanence \(B_r\), the thickness \(d\), the wavenumber \(k = 2\pi/\lambda\) of the wavelength \(\lambda\) and the distance \(y\) from the array. With \(M\) blocks per wavelength, the field is lower by the factor \(\sin(\pi/M)/(\pi/M)\), and a harmonic of order \(M - 1\) reaches the weak side. Mallinson described this one-sided flux, and Halbach used such arrays for the magnets of particle accelerators.

The example solves one wavelength of an endless array as a periodic cell, with four and with eight blocks per wavelength. The field above and below the array must agree with these formulas.

import math

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

The cell

The array has a wavelength of \(\lambda = 40\,\text{mm}\) and is \(d = 10\,\text{mm}\) thick, a quarter of the wavelength, as with square blocks and four blocks per wavelength. Block \(j\) is magnetized at the angle \(360^\circ j/M\). The cell glues the \(M\) blocks to air \(1.5\,\lambda\) deep above and below the array, where the field has faded to almost nothing. Its left and right sides are named left and right and identified, so that the mesh repeats from cell to cell, and the far edges above and below are named far (Figure 1):

lam, d, Br = 40.0, 10.0, 1.3  # the wavelength and the thickness in mm, the remanence in T
k = 2 * math.pi / lam
depth = 1.5 * lam


def cell(M):
    """one wavelength of the array with M blocks, the air above and below, its sides identified"""
    width = lam / M
    blocks = [px.QuickSketch((j * width, 0)).Rectangle(width, d).Face() for j in range(M)]
    for j, block in enumerate(blocks):
        block.name = f"block{j}"
    air = [px.QuickSketch((0, d)).Rectangle(lam, depth).Face(), px.QuickSketch((0, -depth)).Rectangle(lam, depth).Face()]
    for part in air:
        part.name = "air"
    shape = px.Glue(blocks + air)
    for edge in shape.edges:
        (x0, y0), (x1, y1) = edge.bounding_box
        if x1 < 1e-9:
            edge.name = "left"
        elif x0 > lam - 1e-9:
            edge.name = "right"
        elif y1 < -depth + 1e-9 or y0 > d + depth - 1e-9:
            edge.name = "far"
    shape.Identify("left", "right")
    return shape


array = cell(4)
print("faces", [face.name for face in array.faces])

faces ['block0', 'block1', 'block2', 'block3', 'air', 'air']

Figure 1: One wavelength of the array, with the direction of magnetization of each block

The field

With the field \(B = \operatorname{curl} A\) of the vector potential \(A\) along \(z\), and blocks with the permeability of air, the magnetostatic equations become \(\int \nabla A \cdot \nabla v = \int (B_{r,x}\,\partial v/\partial y - B_{r,y}\,\partial v/\partial x)\) for every test function \(v\), with the remanence \((B_{r,x}, B_{r,y})\) of each block. NGSolve’s Periodic space joins the identified sides, \(A\) vanishes on the far edges, and \(B = (\partial A/\partial y, -\partial A/\partial x)\):

def field(M, order=4):
    """the mesh, the vector potential A and the flux density B of the array with M blocks per wavelength"""
    mesh = ngs.Mesh(px.OCCGeometry(cell(M)).GenerateMesh(maxh=lam / 16))
    mesh.Curve(order)
    fes = ngs.Periodic(ngs.H1(mesh, order=order, dirichlet="far"))
    u, v = fes.TnT()
    turn = {f"block{j}": 2 * math.pi * j / M for j in range(M)}
    bx = mesh.MaterialCF({name: Br * math.cos(a) for name, a in turn.items()}, default=0)
    by = mesh.MaterialCF({name: Br * math.sin(a) for name, a in turn.items()}, default=0)
    a = ngs.BilinearForm(ngs.grad(u) * ngs.grad(v) * ngs.dx).Assemble()
    f = ngs.LinearForm((bx * ngs.grad(v)[1] - by * ngs.grad(v)[0]) * ngs.dx).Assemble()
    A = ngs.GridFunction(fes)
    A.vec.data = a.mat.Inverse(fes.FreeDofs(), inverse="sparsecholesky") * f.vec
    return mesh, A, ngs.CF((ngs.grad(A)[1], -ngs.grad(A)[0]))


mesh, A, B = field(4)

In 2D, the field lines are the lines of constant \(A\), and drawn at equal steps of \(A\), they lie the closer together the stronger the field is. They close above the array, and hardly any reach below it:

Against the formulas

On the strong side, the harmonic of order \(M + 1\) adds a ripple that fades as \(e^{-(M+1)ky}\), so a little away from the array, \(|B|\) is the same all along it. On the weak side, only the harmonic of order \(n = M - 1\) is left, with the magnitude \(B_r\,(1 - e^{-nkd})\,\lvert\sin(n\pi/M)/(n\pi/M)\rvert\,e^{-nky}\). The smallest and the largest \(|B|\) along the array at the distances \(\lambda/8\), \(\lambda/4\) and \(\lambda/2\) above and below it are compared with the formulas:

for M in (4, 8):
    mesh, A, B = field(M)
    strong = Br * (1 - math.exp(-k * d)) * math.sin(math.pi / M) / (math.pi / M)
    n = M - 1
    weak = Br * (1 - math.exp(-n * k * d)) * abs(math.sin(n * math.pi / M) / (n * math.pi / M))
    for y in (lam / 8, lam / 4, lam / 2):
        above = [ngs.Norm(B)(mesh(lam * i / 40, d + y)) for i in range(40)]
        below = [ngs.Norm(B)(mesh(lam * i / 40, -y)) for i in range(40)]
        print("M", M, "| y", y, "| above", round(min(above), 5), "to", round(max(above), 5), "formula",
              round(strong * math.exp(-k * y), 5), "| below", round(min(below), 6), "to", round(max(below), 6),
              "formula", round(weak * math.exp(-n * k * y), 6))

M 4 | y 5.0 | above 0.41819 to 0.42752 formula 0.4227 | below 0.035971 to 0.037369 formula 0.036645
M 4 | y 10.0 | above 0.19264 to 0.19282 formula 0.19273 | below 0.003468 to 0.003479 formula 0.003473
M 4 | y 20.0 | above 0.04006 to 0.04006 formula 0.04006 | below 3.1e-05 to 3.1e-05 formula 3.1e-05
M 8 | y 5.0 | above 0.45743 to 0.45766 formula 0.45753 | below 0.000693 to 0.000774 formula 0.000741
M 8 | y 10.0 | above 0.2086 to 0.20861 formula 0.20861 | below 3e-06 to 3e-06 formula 3e-06
M 8 | y 20.0 | above 0.04336 to 0.04337 formula 0.04336 | below 0.0 to 0.0 formula 0.0

The table shows two things. First, the weak side is weak indeed. A quarter wavelength from the array with four blocks, the field is 0.0035 T below it and 0.193 T above it, and with eight blocks the field below falls to a few microtesla. Second, more blocks gain little on the strong side. Eight blocks per wavelength give 0.209 T instead of 0.193 T, 8 % more, since the factor \(\sin(\pi/M)/(\pi/M)\) only rises from 0.900 to 0.974. Most arrays therefore use four blocks per wavelength.

See Periodic identifications for Identify and periodic meshes.