Magnets, electric motors, and transformers are areas where problems involving magnetostatics can be found. The "statics" implies that the time rate of change is slow, so we start with Maxwell's equations for steady cases,

$$\nabla \times H=J$$

$$\nabla \cdot B=0$$

and the relationship

$$B=\mu H$$

where **B** is the *magnetic
flux density,* **H** is the *magnetic
field intensity*, **J** is
the *current density*, and *µ* is
the material's *magnetic permeability*.

Since $$\nabla \cdot B=0$$, there exists a *magnetic
vector potential* **A** such
that

$$B=\nabla \times A$$

and

$$\nabla \times \left(\frac{1}{\mu}\nabla \times A\right)=J$$

The plane case assumes that the current flows are parallel to
the *z*-axis, so only the *z* component
of **A** is present,

$$A=(0,0,A),\text{\hspace{1em}}J=(0,0,J)$$

You can impose the common gauge assumption (Lorenz gauge or
Coulomb gauge, see Wikipedia^{®} [2])

$$\nabla \xb7A=0,$$

and then the equation for **A** in
terms of **J** can be simplified to the
scalar elliptic PDE

$$-\nabla \text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(\frac{1}{\mu}\nabla A\right)=J,$$

where *J* = *J*(*x*,*y*).

For the 2-D case, we can compute the magnetic flux density **B** as

$$B=\left(\frac{\partial A}{\partial y},-\frac{\partial A}{\partial x},0\right)$$

and the magnetic field **H**, in
turn, is given by

$$H=\frac{1}{\mu}B$$

The interface condition across subdomain borders between regions
of different material properties is that **H** x **n** be continuous. This implies the continuity
of

$$\frac{1}{\mu}\frac{\partial A}{\partial n}$$

and does not require special treatment since we are using the variational formulation of the PDE problem.

In ferromagnetic materials, *µ* is usually
dependent on the field strength |*B*| = |∇*A*|,
so the nonlinear solver is needed.

The Dirichlet boundary condition specifies the value of the
magnetostatic potential *A* on the boundary. The
Neumann condition specifies the value of the normal component of

$$n\cdot \left(\frac{1}{\mu}\nabla A\right)$$

on the boundary. This is equivalent to specifying the tangential
value of the magnetic field *H* on the boundary.

Visualization of the magnetostatic potential *A*,
the magnetic field **H**, and the magnetic
flux density **B** is available. **B** and **H** can
be plotted as vector fields.

[1] Popovic, Branko D., *Introductory Engineering
Electromagnetics*, Addison-Wesley, Reading, MA, 1971.

[2] Wikipedia entries on Gauge fixing.

As an example of a problem in magnetostatics, consider determining the static magnetic field due to the stator windings in a two-pole electric motor. The motor is considered to be long, and when end effects are neglected, a 2-D computational model suffices.

The domain consists of three regions:

Two ferromagnetic pieces, the stator and the rotor

The air gap between the stator and the rotor

The armature coil carrying the DC current

The magnetic permeability *µ* is 1 in
the air and in the coil. In the stator and the rotor, *µ* is
defined by

$$\mu =\frac{{\mu}_{\mathrm{max}}}{1+c{\Vert \nabla A\Vert}^{2}}+{\mu}_{\mathrm{min}}.$$

*µ*_{max} = 5000, *µ*_{min} =
200, and *c* = 0.05 are values that could represent
transformer steel.

The current density *J* is 0 everywhere except
in the coil, where it is 1.

The geometry of the problem makes the magnetic vector potential *A* symmetric
with respect to *y* and antisymmetric with respect
to *x*, so you can limit the domain to *x* ≥
0,*y* ≥ 0 with the Neumann boundary condition

$$n\text{\hspace{0.17em}}\cdot \text{\hspace{0.17em}}\left(\frac{1}{\mu}\nabla A\right)=0$$

on the *x*-axis and the Dirichlet boundary
condition *A* = 0 on the *y*-axis.
The field outside the motor is neglected leading to the Dirichlet
boundary condition *A* = 0 on the exterior boundary.

The geometry is complex, involving five circular arcs and two
rectangles. Using the PDE app, set the *x*-axis limits
to [-1.5 1.5] and the *y*-axis limits to ```
[-1
1]
```

. Set the application mode to **Magnetostatics**,
and use a grid spacing of 0.1. The model is a union of circles and
rectangles; the reduction to the first quadrant is achieved by intersection
with a square. Using the "snap-to-grid" feature, you
can draw the geometry using the mouse, or you can draw it by entering
the following commands:

pdecirc(0,0,1,'C1') pdecirc(0,0,0.8,'C2') pdecirc(0,0,0.6,'C3') pdecirc(0,0,0.5,'C4') pdecirc(0,0,0.4,'C5') pderect([-0.2 0.2 0.2 0.9],'R1') pderect([-0.1 0.1 0.2 0.9],'R2') pderect([0 1 0 1],'SQ1')

You should get a CSG model similar to the one in the following plot.

Enter the following set formula to reduce the model to the first quadrant:

(C1+C2+C3+C4+C5+R1+R2)*SQ1

In boundary mode you need to remove a number of subdomain borders.
Using **Shift**+click, select borders and remove them
using the **Remove Subdomain Border** option
from the **Boundary** menu until the geometry
consists of four subdomains: the stator, the rotor, the coil, and
the air gap. In the following plot, the rotor is subdomain 1, the
stator is subdomain 2, the air gap is subdomain 3, and the coil is
subdomain 4. The numbering of your subdomains may be different.

Before moving to the PDE mode, select the boundaries along the *x*-axis
and set the boundary condition to a Neumann condition with *g* =
0 and *q* = 0. In the PDE mode, turn on the labels
by selecting the **Show Subdomain Labels** option
from the **PDE** menu. Double-click each subdomain
to define the PDE parameters:

In the coil both

*µ*and*J*are 1, so the default values do not need to be changed.In the stator and the rotor

*µ*is nonlinear and defined by the preceding equation. Enter*µ*as5000./(1+0.05*(ux.^2+uy.^2))+200

`ux.^2+uy.^2`

is equal to |∇*A*|^{2}.*J*is 0 (no current).In the air gap

*µ*is 1, and*J*is 0.

Initialize the mesh, and continue by opening the Solve Parameters
dialog box by selecting **Parameters** from
the **Solve** menu. Since this is a nonlinear
problem, the nonlinear solver must be invoked by checking the **Use
nonlinear solver**. If you want, you can adjust the tolerance
parameter. The adaptive solver can be used together with the nonlinear
solver. Solve the PDE and plot the magnetic flux density *B* using
arrows and the equipotential lines of the magnetostatic potential *A* using
a contour plot. The plot clearly shows, as expected, that the magnetic
flux is parallel to the equipotential lines of the magnetostatic potential.

**Equipotential Lines and Magnetic Flux in a Two-Pole Motor**

Was this topic helpful?