Open-boundary FEA without a giant air box

Every open-boundary magnetostatic model — a magnet in free space, a solenoid actuator, an unshielded machine, a pair of coils — ends at some artificial outer boundary. What you put on that boundary decides how much of your mesh is device and how much is air you never wanted to solve.

The two boundaries everyone starts with are both wrong

A Dirichlet condition, A = 0 on the outer edge, forces every flux line to close inside the box. It under-estimates the reluctance of the exterior path, and the error decays only as the box grows. A Neumann condition, ∂A/∂n = 0, forces flux lines to cross the boundary at right angles — equally wrong, with the opposite sign of error.

Neither is a mistake exactly; both converge to the right answer as the box grows. The problem is the price. A “good enough” Dirichlet box is typically 10 to 20 times the device radius, which means the overwhelming majority of your elements are solving for air you have no interest in — and every one of them costs you in assembly, in factorisation, and in a transient run at every single time step.

What the exterior field actually looks like

Outside every source, A satisfies Laplace's equation, and the general exterior solution that decays at infinity is a series of harmonics:

A(r, θ) = a₀ ln r + Σₙ≥₁ r⁻ⁿ (aₙ cos nθ + bₙ sin nθ)

Take one harmonic n. On the circle r = R it satisfies ∂A/∂r = −(n/r)·A exactly, which rearranges into a homogeneous Robin condition:

∂A/∂n + (n/R)·A = 0        on r = R

That is the far-field, asymptotic or “balloon” boundary condition. It is not an approximation of the exterior field — it is exact for a pure r⁻ⁿ exterior harmonic. Its residual error is the amplitude, at r = R, of the harmonics you did not name, each of which is already decaying as (r_device/R)ᵐ.

Why it is nearly free to implement, and free to use

Substituted into the weak form, the condition contributes a pure stiffness term and nothing at all to the right-hand side:

Kᵢⱼ  +=  (n / (μ₀ R)) · ∮ Nᵢ Nⱼ dl

That is one symmetric edge-mass block per boundary segment. No extra mesh, no extra unknowns, no new element type. The added block is symmetric and positive semi-definite, so the matrix stays SPD and every solver in the code — conjugate gradient, sparse LU, the Newton tangent — is unchanged in character. Because the contribution to the load vector is exactly zero, the Newton residual picks the term up with the right sign on its own.

The alternatives are considerably more invasive. A Kelvin transformation (“ballooning” proper) meshes a second, inverted disc glued one-to-one to the first: essentially exact, but it doubles the air mesh and needs matched node counts on two circles. Infinite elements need new shape functions in every assembly kernel. The Robin condition buys most of the accuracy for a fraction of the machinery.

What it measures out at

Planar dipole, against a converged Dirichlet reference at R = 50 r_device. At R = 3 r_device the far-field boundary reproduces that reference to 0.09 % in A; a Dirichlet box at the same radius is off by 12.5 %. The axisymmetric variant matches to 0.03 % on flux linkage.

The axisymmetric case shows the shape of the effect most clearly. Maximum error in u = r·Aϕ over eight off-axis air probes, current-loop model, against a Dirichlet reference at R = 50 r_device:

R / r_device	Dirichlet	Far-field
1.5	40.0 %	0.90 %
2	17.1 %	0.13 %
3	5.05 %	0.06 %
5	1.07 %	0.05 %

The Dirichlet column is converging, slowly, exactly as theory says it must. The far-field column has reached its floor by R = 2. In practice a boundary at two to three device radii with the right decay order replaces a Dirichlet box at twenty — and in 2D that is roughly two orders of magnitude of air elements you no longer pay for.

The two things you have to get right

1. The decay order n

This is the one genuinely physical input: which exterior harmonic dominates at your boundary.

Model	Leading harmonic	n
Magnet, coil pair, actuator, conductor pair	dipole	1 (default)
Quadrupole, e.g. a 4-pole focusing magnet	quadrupole	2
Full cross-section of a p-pole-pair machine	p-th	p
Symmetry sector with periodic boundaries	p-th	p

n = 1 is correct, or conservative, for almost everything — which is why it is the default. Getting it wrong does not blow the model up; it leaves you with a boundary that behaves like a mildly improved Dirichlet.

2. The net current has to be zero

If the modelled domain carries a non-zero net current then a₀ ≠ 0 and A ~ ln r — a field that no r⁻ⁿ condition can represent, and that “A = 0 at infinity” does not even admit. The remedy is not a different boundary condition; it is to put the return conductor in the model. Nabla's checklist warns when a far-field boundary coexists with a net source current rather than letting it pass quietly.

Two details that bite in axisymmetric models

In axi mode the unknown is the flux function u = r·Aϕ, the exterior expansion is in spherical coordinates, and R is the spherical radius ρ = hypot(r, z − z_c) measured from the field centre — not the cylindrical radius. So an axisymmetric far-field boundary wants to be drawn as a hemispherical arc centred on the field centre and terminating on the axis, not as a cylinder.

The axi weak form also carries a 1/r weight, so the boundary integral is 1/r-weighted — the opposite r-power from a thermal convection (Robin) boundary on the same mesh. The table above is the measurement that settled which one is right: built with +r weighting instead, the same condition merely behaves like a slightly improved Dirichlet (19.8 % at R = 1.5 against 0.90 %; 0.57 % at R = 5 against 0.05 %), which is precisely what an under-weighted Robin coefficient should do. Two solvers, identical except for that one factor, is a cheaper way to settle the question than an argument about it.

In Nabla

The boundary type is called Far-field (balloon), and it takes a decay order, a field centre, and a radius R where 0 means auto. Auto is the recommended setting: R is then derived per segment from the mesh, which removes the most common user error — a hand-entered R inconsistent with the drawn geometry — and degrades gracefully on a boundary that is not a perfect circle. A “centre on boundary centroid” button fills the centre in from the selected shapes.

It composes with everything else, because it is only an addition to K: Newton iteration, Crank–Nicolson time integration, time-harmonic solves, monolithic circuit coupling, periodic and anti-periodic sectors, the sliding band, and second-order elements. Far-field boundaries are part of Nabla Core.

Related: second-order (P2) elements, several motion zones in one model, and the boundary conditions chapter of the manual.