**By R. Belmans**

## Enhanced accuracy of finite element field quantities

It will be focused on the practical application of the static electromagnetic field solution of Laplace’s equation in a local post-process to increase the accuracy of an existing solution obtained by the standard finite element method using first order elements. Advantages and drawbacks are discussed.

**The Laplace equation****in two dimensions, expressed in polar co-ordinates****is:**

**Assuming linearity and uniformity, and applying a Fourier series to eq.(5.317), yields die harmonic function:**

with its coefficients:

The procedure for solving eq.(5.318) describes the solution of a Dirichlet problem on a circle with given boundary values at its circumference. The coefficientscan be calculated using known potentialsat the circumference of a circle with radius R.

Now a finite number of N equi-angularly arranged points are applied onto the circumference of the circle.

With N boundary potential valuesknown on the circumference and according to the properties of harmonic functions the first term in eq. (5.318) can be written by:

**The Fourier coefficients are rewritten as follows:**

With the Fourier series (5.318) and their coefficients eq.(5.322) the potential in the centre of a circle can be computed knowing only the boundary potential values on the circumference of the circle.

**Using this approach inside a finite element solution,** the value of the potential of a field point now depends on the solution in several finite elements. Thus, local numeric errors in single elements have a relatively small influence on the solution in the considered field point. Applying (5.318) derivatives at the centre of the circle can be calculated in a closed analytical form, avoiding numerical differentiation.

**The idea is to adapt the described process,** of solving a Dirichlet problem on a circular surface, to determine the vector potential in a point Pi of a discretised finite element domain (Fig. 5.96). Ri is the radius of the considered circular surface and the dots at the circumference indicate the points of known vector potential values computed beforehand. These points do not have to be nodes of the actual finite element mesh. This feature makes the technique very advantageous to automatic and adaptive meshing schemes in which the user can not guarantee the control of the mesh and especially its symmetry,

To obtain the potential distribution at a given contour inside a finite element domain, multiple circles have to be evaluated. Overlapping circles guarantee a continuous solution in the considered region after the post-process.

**The numerical shape of (5.318), (5.321) and (5.322)** enables an easy implementation of the procedure in a finite element program package. The derivatives in the centre of the circle are represented by the Fourier coefficients. Thus, no additional computational effort is necessary to compute the flux density in the centre.

**Fig. 5.96. Multiple circles to determine the vector potential on a contour.**

**The local solution of the** Laplace equation inside an air gap of an electromagnetic device, using a Fourier series approximation for the vector potential, results in a significant increase in accuracy of the derived field quantities. To compare the results obtained by the local field evaluation to the conventionally obtained field quantities of first order elements, Fig. 5.97 shows the computed magnetic flux density derived by using the Laplace approach. For this application of the local Dirichlet problem, 24 potential boundary values on the circumference of the circle were used.

**Applying such field quantities** to the Maxwell stress tensor to compute the local forces acting on bodies inside an electromagnetic field yields values of higher accuracy. To obtain the torque of an electrical machine, the local force values are integrated along a contour in the airgap-

**Fig. 5.97. a) Computed vector potential inside the circular FEM domain and b) the resulting flux density Bx derived by applying the local post-process.**

**Another approach to compute** the torque more accurately, uses the values of the magnetic vector potential on two concentric circles with radiias boundary conditions (Fig. 5.98). Local field values on the circular contour C with radiusare calculated.

**Fig. 5.98. Local Dirichlet problem for a cylindrical air gap.**

**If the inner radius****is taken as a reference, the general solution of Laplace’s equation is:**

**The coefficients****are independently determined for each circular harmonic. A fast Fourier transformation (FFT) algorithm is used to express the magnetic vector potential at the boundaries as a series of such circular harmonics:**

**Once the magnetic vector potential at the contour C is known, the normal and tangential component of the magnetic flux density can be determined:**

The tangential force componentresults in the torque T of the device, ft can be shown (Salon Mertens et al. 7S) that the value of the torque is given by

being independent of the radius r of contour C. It is not necessary to calculate the normal and tangential component of the magnetic flux density on the contour resulting in a faster algorithm, when the overall torque is aimed at. The proposed method can easily be extended to time-harmonic problems. If all values are rms-values the torque is obtained by adding the torque calculated using the real- and the imaginary-component of the solution.

This method has its advantage for electrical machine analysis as it is suited for small air gaps. By using eq.(5.329) the torque is evaluated directly, without the explicit calculation of the flux densities.

**Fig. 5.99. Equipotential plot of the real component solution of a 400 kW induction motor.**

The performance of this method is compared with the classical Maxwell stress tensor method using a model of a 400 kW induction machine for tests.

**Table 5.10. Data of the 400 kW induction machine.**

To ensure accurate results, a good trade-off between mesh refinement and using the enhanced post-processing methods is necessary. A relatively coarse discretisation in the air gap of the induction machine was chosen (Fig. 5.100).

**With such a coarse discretisation,** the Laplace-based method is less sensitive to the actual choice of the contour inside the air gap than the classical method. The air gap spans a region between an inner radius of 0.186 m up to an outer radius of0.1875 m. Fig. 5.101 shows the variation of the calculated torque. Contours with different radii are chosen. For the Laplace-based method, the inner and outer radii are varied simultaneously. Therefore, the value of the torque varies symmetrically towards the middle of the air gap.

**Fig. 5.100. Detail of the discretisation in the air gap.**

**The variation of the calculated torques** using the enhanced method is much smaller when compared to the classical method. It must be stated, however, that an appropriate mesh refinement scheme would lead to better results even for the classical torque computation.

**Fig. 5.101. Variation of the torque calculated along different contours inside the air gap.**

**The rate of convergence of the relative** error indicates that the required number of sample points for the enhanced method can be chosen to be substantially fewer than for the classical method. As shown in Fig. 5.102, the smallest relative error is computed with 16384 points using the classical Maxwell stress method, instead of 2048 points by using the Laplace-based approach.

**Fig. 5.102. Rate of convergence of the relative error.**

**The same basic idea, as used in the ‘circle’ approach**, yields the local solution for the three-dimensional field. The local field problem is now defined by the known potential values equally distributed along the surface of a sphere assumed to be the boundary potential values of the local field problem. According to the co-ordinate transformation:

a spherical co-ordinate system is applied (Fig. 5.103).

**Fig. 5.103. a) Sphere with b) co-ordinate system.**

Using the Laplace equation with the co-ordinate transformation yields:

**Applying the theorem of the separation of the variables ****to (5.330), a general form of the functions ****depending on the potential A can be written. Every solution of the Laplace equation, being finite for all****is a solution of:**

associated Legendre polynomial of the first kind. To simplify the notations, the surface harmonics

are introduced. Assuming (5.331) to be a linear form, the potential in the origin isfinite. The constantsare linear combinations ofThe summation

is a solution of (5.330). Here, the magnetic scalar potential A is completely determined by the constantsThe aim is to calculate the magnetic flux density at a point using known scalar potential values in its vicinity. Consequently a spherical volume with known boundary potentials at its surface around this field point is chosen to determine the field. The known boundary potentials result, as in the two-dimensional case, from a previously performed FEM computation and determine all constants in (5.333). To calculate the magnetic field quantities at the centre of the spherical volume, the Laplace equation has to be solved locally and spherically around this field point with radius The boundary potential values are available only as single values at the surface of the sphere. To distribute them equally along this surface, the spherical co-ordinatesare divided into /andinto K equal angles

**In order to satisfy (5.333) accurately**, the numbers of J and K must be sufficiently large. On the other hand, large numbers increase the computational expenses rapidly. With respect to the computation time and accuracy a good compromise has to be found. Practical values for J and K are given in the following section.

**Assuming****to be the co-ordinates in the local system with the interesting field point at the centre of the sphere, the coefficients**** and****can be determined by a Legendre decomposition using the boundary potential values:**

**Retaining the local co-ordinate system in****and the magnetic flux density in the original global co-ordinate system, and calculating the derivatives at the origin of the local co-ordinate system (Fig. 5.103) using****in (5.333), yields:**

Analogous to (5.336) the derivatives, inare found by taking and inby takingin (5.333). With respect to (5.336), using (5.335) and with the Legendre terms:

the components of the flux density at the centre of a sphere are explicitly rewritten by:

**Using this local field approach (5.338)** by arranging multiple overlapping spheres at an arbitrary surface or contour (Fig. 5.104), it is possible to obtain the required local field quantities at this surface with the same accuracy as the previously by FEM computed potential values. The tetrahedron shown in Fig. 5.104 represents a part of the three-dimensional mesh of the FEM domain.

**Fig. 5.104. Arrangement of multiple overlapping spheres to obtain the local field values on an arbitrary contour/surface across the centre points of the spheres inside a three-dimensional FEM domain.**

**Fig. 5.105. Flux density distribution****on the front surface of****(see Fig. 5.106), a) computed by the classical direct derivation of the potential and b) using the proposed post-processor method.**

**Fig. 5.106. Three-dimensional FEM model of the test example.**

**From Fig. 5.**105 the difference between the direct evaluation of the potential and the new post-process operator is shown. Here Bz is computed for a test example (Fig. 5.106) at the front surface offacing the permanent magnet cube. It is obvious that linear shape functions, approximating the scalar potential, result in a piecewise constant flux density distribution (Fig. 5.105a). Computed forces starting from this type of solution are unreliable. The local values of plotted in Fig. 5.105b show the expected continuous distribution computed using the new post-processor method.

**Fig. 5.107. Comparison of the convergence behaviour of the FEM potential solution with both the direct derivation and the derivative-free approach.**

**In Fig. 5.107 the quadratic convergence,** referred to the characteristic length h of a finite element, of the FEM potential solution and the rate of convergence of the force computations using both the classical and the new post-processing approach, is plotted versus the number of tetrahedron elements. The same statement can be made for the two-dimensional approach. The triangles in Fig. 5.107 indicate the theoretical gradient of convergence eq.(5.307). The refinement of the three-dimensional discretisation is performed in such a way that the elements are of the same shape in every FEM model to obtain a regularly distributed mesh for all cases. To compute the total force, the Maxwell stress tensor is used integrating the force density calculated in points equidistantly distributed by the density D on all six sides of T. For the classical approach a density D=40 is chosen and in the case of the new method, D is set to 7. The sphere parameters are J-K~\5. The integration surface of the force computations is located in such a way that no plane of r cuts through the nodes of the FEM mesh. If nodes coincide with the points of the force computation using the classical post-processor approach, this would result in a larger error due to the troublesome definition of normal and tangential field components in a node of an element. The gradient-triangles in Fig. 5.107 indicate the theoretical rate of convergence for the quadratic and the linear convergence case. It can be seen as theoretically expected that the relative error in an energy norm of the FEM potential solution converges quadratically, referred to the specific diameter h of the elements eq.(5.307), by increasing the number of first order tetrahedron elements. Due to the analytically described potential function inside the local field volumes, the resulting overall force using this approach is of the same order of convergence. Therefore, no loss of accuracy of the derived field quantities occurs. The convergence of the total forces, computed by the classical approach, indicates the expected linear behaviour. The accuracy of the computed values is influenced by the numerically-obtained derivatives. This shows that the results obtained by the classical method are inherently inaccurate when compared to the accuracy of the potential solution.

**The use of the proposed approach** to enhance the accuracy of computed field quantities starting from an existing potential solution demands an additional step during the post-processing of the FEM analysis (Fig. 5.108).

**Fig. 5.108. Additional step during post-processing to enhance the accuracy of derived field quantities.**

**Having obtained a FEM potential solution**, the user only has to define the surface of integrationon which the field quantities or forces have to be calculated. Defining an arbitrary contour allows the computation of field quantities or forces along it as well. For each plane or contour the density D, the sphere parameter J, K and the radius R have to be set. The sphere parameter are problem-dependent and related to the geometry of the device, i.e. the air gap width. The planes or contours should be centred in the air gap. A suitable value for the diameter of the single spheres is about 90-95 % of the air gap width to have as many tetrahedron finite elements inside the sphere as possible. Including only one finite element in the sphere results in no enhancements in accuracy of the derived quantity. To ensure a continuous field solution, the density D should be chosen in such a way that the spheres overlap (Fig. 5.104). For the distance between two points on the surface of integration, it is suitable to choose the radius of the sphere.

**Fig. 5.109. Spheres with different numbers for the parameter J, K.**

**To define the number and position** of boundary potential values distributed on the surface of each sphere, the parameters J and K have to be chosen. To ensure uniformly distributed boundary values J is set equal to K. In accordance with the results of Fig. 5.105 and other test calculations a numberis sufficient to meet the ratio between computational costs and accuracy. Fig. 5.109 illustrates by different J=K the position and number of boundary potentials to approximate the local field inside a sphere.