Numerical implementation of the FEM (Electrical Machine) Part 1


In this section, the theoretically derived abstract knowledge of the finite element method will be practically applied to magnetostatic field problems. Standard linear triangular finite elements will be used. The system of equations will be derived using the energy minimum functional.

Starting with the numerical solution of the Laplace equation in two dimensions, then introducing impressed currents and permanent magnet material, the difficulties programming a FEM code will be introduced and treated. Reading this section should enable the reader to obtain an understanding of the practical realization of the method in a computer program.

Laplace’s equation

For simplicity in this example of the implementation of the FEM, a two-dimensional approach is considered. A Cartesian co-ordinate system is assumed. In a two-dimensional field the electrical field strength consists of one component in the z-direction only and the magnetic fields are in the xy-plane.


The vector potential A and the current density J have a component only in the z-direction. For reasons of simplicity the index z, to indicate the direction, will not be written in the following.


The static magnetic field in terms of the vcctor potential is given by the A-formulation of Laplace’s equation:


At this moment air is considered as material. Using the energy minimum functional


For a complete problem definition a magnetic flux has to be imposed by applying the appropriate Dirichlet conditions at the boundaries of the field domain of interest.

Linear basis function

The problem region is discretised by triangular node elements. The vector potential is approximated by linear shape functions. By connecting all triangular elements at their nodes and forcing the node potential to be equal, the magnetic vector potential becomes continuous over the defined field region. Linear shape functions are applied by using:


The coefficients a, b and c are found from the values of the magnetic vector potentialtmp10867_thumbat the three nodes of an element (Fig. 5.42).






Triangular finite element.

Fig. 5.42. Triangular finite element.

The magnetic vector potential within an element is:


Using the shape functions, the magnetic vector potential is approximated by:




The shape functions _tmp10875_thumb_ are found by cyclic permutation of the indices in (5.101) and (5.102). The areatmp10876_thumbof an clement is:


It is easy to verify that the three shape functions are one at one node and zero at both other nodes.


Functional within an element

The gradient of the magnetic vcctor potential in terms of the shape functions can be calculated by:


Substituting (5.105) into the functional within an element gives


The elements of the 3 x 3 element matrixtmp10885_thumbare


In a matrix-vector notation, the functional within an element is:


Assembling all elements

The overall functional is found as the sum of the functionals within the elements.


Summation is practically done element by element. The element matrix Kw is calculated for each finite element and added to the already existing coefficient matrix K. For the two finite elements in Fig. 5.43 this yields:


Two elements combined with each other have four nodes, i.e. the coefficient matrix K becomes a 4 x 4 matrix.

 Joining two finite elements.

Fig. 5.43. Joining two finite elements.

The system of equation

The system of linear equations is found by forcing the partial derivatives with respect to all unknowns to be zero. This is in fact an energy minimisation.


Applying Dirichlet conditions at boundaries is fixing the magnetic vector potential at some prescribed values. Two groups of nodes can be distinguished: nodes with free potentials that are to vary (index f) and nodes with prescribed potential values (index p). The node numbering is such that all nodes with potentials that arc free to vary are numbered first.


Applying the boundary conditions, the system of linear equations becomes:


The generated coefficient matrix owns some interesting properties:

• Sparse,tmp10896_thumbis different from zero only if nodetmp10897_thumbandtmp10898_thumbare connected by an element.

• The matrix is symmetrictmp10899_thumbdiagonal dominant and positive definite.

The assembly of the matrix is straightforward and is done element by element (Fig. 5.44).Assembling of the coefficient matrix and the right hand side.

Fig. 5.44. Assembling of the coefficient matrix and the right hand side.

Material properties

Laplace’s equation in two dimensions is given by


With the reluctivity, respectively permeability and their relative values, index


characterizing the ferromagnetic properties of the material used. Due to numerical reasons the relative reluctivitytmp10907_thumbis used:



The functional is:


and is valid for magnetostatic problems with different materials and where only Dirichlet and Neumann conditions are applied.

Functional within an element

The functional within an element becomes:



The elements of the element matrix


Poisson’s equation

Poisson’s equation in two dimensions is given by:


or due to numerical reasons


The functional is


and is valid for magnetostatic problems with different materials and given current densities. Dirichlet and Neumann conditions can be imposed at the boundaries.

Functional within an element

The functional within an element becomes


The evaluation of the second term gives:


In matrix-vector notation, the functional within an element is


The elements of the source vectortmp10921_thumbare found as


Solution of Poisson’s equation

The system of linear equation is found by forcing the partial derivatives with respect to all unknowns to be zero.