5.2 Three-Dimensional Vector Finite Element Method

The vector finite element analysis will be demonstrated using the quasi magnetostatic case as described in Subsection 5.1.3. The Galerkin method is applied to (5.15) and the weak formulation for the vector differential operator, as discussed in Section 3.3, is used. This gives the following expression

$$\begin{split}& \quad \int_{\mathcal{V}}\left(\vec{\nabla}\tim... ...thcal{V}}\mu\vec{N}_i\cdot\vec{H} \mathrm{d}V = 0. \end{split}$$ (5.17)

Adding an arbitrary gradient field (i. e. $\vec{\nabla}\psi$ ) to the magnetic field $\vec{H}$ does not alter (5.17), since the rotor operator of a gradient field is zero ( $\vec{\nabla}\times\vec{\nabla}\psi{=}0$ ). The solution will remain unchanged like

$$\begin{split}& \int_{\mathcal{V}}\left(\vec{\nabla}\tim... ...thcal{V}}\mu\vec{N}_i\cdot\vec{H} \mathrm{d}V = 0. \end{split}$$ (5.18)

Applying the substitution

$$\vec{H} + \vec{\nabla}\psi = \vec{H}_1$$ (5.19)

in (5.18) leads to the following equation

$$\begin{split}&\quad \int_{\mathcal{V}}\left(\vec{\nabla}\time... ...dot(\vec{H}_1 - \vec{\nabla}\psi) \mathrm{d}V = 0. \end{split}$$ (5.20)

In the literature the vector potential $\vec{H}_1$ of the current density $\vec{J }$ is often denoted as T and the auxiliary scalar field $\psi$ is denoted as $\Omega$ or $\phi$ , respectively. This gives the widely used names of the numerical procedure for handling the quasi-static method -- T-$\Omega$ or T-$\phi$ method [70,71,72,73,74,75,76,77,78,79]. Using this technique the finite element method can be combined with the boundary integral method to diminish the number of unknowns taking into account the unbounded regions [80]. Of coarse, the dominant magnetic field model can be also considered from the equations for the electrodynamic potentials from Subsection 5.1.2 to obtain expressions for the magnetic vector potential $\vec{A}$ and the electric scalar potential $\varphi$ for the quasi-magnetostatic case [81].

In this chapter the finite element analysis with vector shape functions for the approximation of $\vec{H}_1$ and with scalar shape functions for the approximation of $\psi$ is comprehensively explained. The same can be applied for the electric field $\vec{E}$ from Subsection 5.1.1 or for the electrodynamic potentials $\vec{A}$ and $\varphi$ from Subsection 5.1.2.

De facto (5.20) corresponds to the partial differential equation

$$\vec{\nabla}\times\left(\frac{1}{\gamma}\vec{\nabla}\times\vec{H}_1\right) + \jmath\omega\mu(\vec{H}_1 - \vec{\nabla}\psi) = \vec{0}.$$ (5.21)

Since there are two unknowns in (5.21) -- the vector field $\vec{H}_1$ and the scalar field $\psi$ , an additional relation between these two fields is required. Equation (5.20) (or the equivalent one (5.21)) is derived from (5.11) and (5.13), but (5.12) has not been used. With (4.7) it leads to

$$\vec{\nabla}\cdot\left[\mu(\vec{H}_1 - \vec{\nabla}\psi)\right] = 0.$$ (5.22)

Equation (5.22) provides the additional relation between $\vec{H}_1$ and $\psi$ . Thus the magnetic field $\vec{H}$ is obtained by solving the boundary value problem, given by the partial differential equation system consisting of (5.21) and (5.22) for $\vec{H}_1$ and $\psi$

$$\vec{H} = \vec{H}_1 - \vec{\nabla}\psi.$$ (5.23)

In ideal dielectric regions (5.21) cannot be used, because $\gamma = 0$ . In such regions the current density is zero and $\vec{H}$ can be expressed as gradient field like

$$\vec{\nabla}\times\vec{H} = \vec{0} \Rightarrow \vec{H} = \vec{\nabla}\psi.$$ (5.24)

Thus for the numerical analysis (5.24) can be used in the dielectric part of the simulation domain and (5.21), (5.22), and (5.23) can be applied in the remaining parts, respectively. However (5.24) is valid only for simply connected regions. For non-contiguous regions specific cutting algorithms have to be addressed [82,83,84,85]. Unfortunately, these algorithms are quite expensive pre-processing steps for complex structures. In this work it is preferred to assume low conductivity of the dielectrics and to use (5.21), (5.22), and (5.23) in the entire simulation domain. The conductivity in the dielectric regions should be sufficiently low with respect to the conductivity of the conducting areas. However, it should not be very low, because the first summand of the right hand side of (5.21) would increase dramatically and would lead to an extremely ill conditioned linear equation system.

The boundary $\partial\mathcal{V}$ of the calculation domain $\mathcal{V}$ is divided into a Dirichlet boundary $\mathcal{A}_{D1}$ for $\vec{H}_1$ and a Neumann boundary $\mathcal{A}_{N1}$ for (5.21) and into a Dirichlet boundary $\mathcal{A}_{D2}$ for $\psi$ and a Neumann boundary $\mathcal{A}_{N2}$ for (5.21).

$$\partial\mathcal{V} = \mathcal{A}_{D1} + \mathcal{A}_{N1} \mathrm{and} \partial\mathcal{V} = \mathcal{A}_{D2} + \mathcal{A}_{N2}.$$

The weighting of the residual of (5.21) and the following weak formulation for the finite element analysis has been already dealt with (5.20). The same must be done for (5.22) as well. For this purpose the residual of (5.22) is weighted by a set of scalar functions $\lambda_i$

$$\int_{\mathcal{V}}\vec{\nabla}\cdot\left[\mu(\vec{H}_1 - \vec{\nabla}\psi)\right]\lambda_i \mathrm{d}V = 0$$ (5.25)

and then the first scalar Green's theorem is applied similarly as shown in Section 3.3 for the scalar differential operator. This leads to the equation

$$\int_{\partial\mathcal{V}}\vec{n}\cdot\left\{\lambda_i\left[\mu(\... ...}\lambda_i\cdot\left[\mu(\vec{H}_1 - \vec{\nabla}\psi)\right] \mathrm{d}V = 0.$$ (5.26)

As usual for the finite element analysis the unknown vector function $\vec{H}_1$ is approximated by a sum of known vector functions multiplied by coefficients

$$\vec{H}_1 \simeq \sum_{j=1}^{m}c_j\vec{N}_j + \vec{v}, \mathrm{with} \vec{v} = \sum_{j=n+1}^{M}c_j\vec{N}_j.$$ (5.27)

The same is done for the scalar function $\psi$

$$\psi \simeq \sum_{j=m+1}^{n}c_j\lambda_j + v, \mathrm{with} v = \sum_{j=M+1}^{N}c_j\lambda_j.$$ (5.28)

The functions $\vec{v}$ and $v$ comply with the Dirichlet boundary conditions for $\vec{H}_1$ and $\psi$ , respectively. The tilde sign written over the approximated field quantity is not used any more. Of coarse it will be kept in mind that the sums (5.27) and (5.28) are approximations. As already mentioned the coefficients $c_j$ in the scalar approximation (5.28) correspond to the scalar field values in the nodes of the discretized domain. It will be shown in Subsection 5.2.2 that the coefficient $c_j$ in the vector approximation (5.27) complies with the tangential component of the field $\vec{H}$ along the global edge $j$ of the discretized domain. Similarly to the nodes, the edges are numbered globally for the entire domain and locally in each element. Again there is a connectivity array which binds the local (element) edge index with the corresponding global one. The coefficients for both, nodes and edges lying at the Dirichlet boundary, can be obtained from the Dirichlet boundary condition and are arranged behind the unknown coefficients in the following way

$$j = \underbrace{\overbrace{1........................................ ...es} \overbrace{..................N}^{Dirichlet nodes}}_{known coef\!ficients}$$

The unknown coefficients associated with the edges are numbered from $1$ to $m$ and the unknown ones associated with the nodes count from $m{+}1$ to $n$ , respectively. The known coefficients for the edges given by the Dirichlet boundary of (5.27) are numbered from $n{+}1$ till $M$ and the known coefficients for the nodes on the Dirichlet boundary of (5.28) are numbered from $M{+}1$ to $N$ . Thus the first $n$ coefficients are unknown and are the solution of the finite element calculation and the remaining coefficients numbered from $n{+}1$ to $N$ are given by the Dirichlet boundaries.

Each global vector function $\vec{N}_j$ corresponds to the global edge $j$ and is constructed from element functions $\vec{N}_j^e$ similarly with the scalar functions $\lambda_j$ . In Subsection 5.2.2 it will be demonstrated that $\vec{N}_j$ has no tangential component along other edges except edge $j$ . Thus the boundary term in (5.20) can be given in a different form as

$$\begin{split}& \quad \int_{\partial\mathcal{V}}\vec{n}\cdot\l... ...\cdot\left(\vec{E}\times\vec{n}\right) \mathrm{d}A \end{split}$$ (5.29)

with

$$\vec{E} = \frac{1}{\gamma}\vec{J} = \frac{1}{\gamma} \vec{\nabla}\times\vec{H}_1.$$

Since the global edge functions $\vec{N}_i$ with $i\in[1;m]$ have no tangential component on the edges lying on the Dirichlet boundary $\mathcal{A}_{D1}$ , they must be perpendicular to $\mathcal{A}_{D1}$ (or parallel to $\vec{n}$ ). As clearly shown by the third member of (5.29) the boundary integral in (5.20) has a contribution only for the Neumann boundary $\mathcal{A}_{N1}$ . Furthermore the following is assumed: The electric field $\vec{E}$ is either perpendicular to $\mathcal{A}_{N1}$ , which means that $\vec{E}\times\vec{n}$ is zero in the last term of (5.29), or the simulation domain is sufficiently large and allows that $\vec{E}\times\vec{n}$ can be neglected on $\mathcal{A}_{N1}$ . Accepting this, the boundary term in (5.20) vanishes

$$\int_{\partial\mathcal{V}}\vec{n}\cdot\left[\vec{N}_i\times \left... ...al\mathcal{V}}\vec{N}_i\cdot\left(\vec{E}\times\vec{n}\right) \mathrm{d}A = 0.$$ (5.30)

The boundary integral of (5.26) can be written in the following way

$$\int_{\partial\mathcal{V}}\vec{n}\cdot\left\{\lambda_i\left[\mu(\... ...mathrm{d}A = \int_{\mathcal{A}_{N2}}\lambda_i \vec{n}\cdot\vec{B} \mathrm{d}A$$ (5.31)

with

$$\mu(\vec{H}_1 - \vec{\nabla}\psi) = \mu\vec{H} = \vec{B}.$$

The global node functions $\lambda_i$ are non-zero only in the neighbor elements of the unknown nodes, which do not belong to the Dirichlet boundary $\mathcal{A}_{D2}$ and are indexed by $i\in[m{+}1;n]$ . On the Dirichlet boundary $\mathcal{A}_{D2}$ the global functions $\lambda_i$ are defined to zero. Thus the boundary integral (5.31) vanishes at the Dirichlet boundary $\mathcal{A}_{D2}$ and the integration domain is restricted to the Neumann boundary $\mathcal{A}_{N2}$ . Additionally it is assumed that the magnetic flux $\vec{B}$ is perpendicular to $\vec{n}$ on the surface $\mathcal{A}_{N2}$ or that it can be neglected on this surface. Thus the boundary integral of (5.26) is also set to zero

$$\int_{\partial\mathcal{V}}\vec{n}\cdot\left\{\lambda_i\left[\mu(\... ...}A = \int_{\partial\mathcal{V}}\lambda_i \vec{n}\cdot\vec{B} \mathrm{d}A = 0.$$ (5.32)

In such case it is often spoken of a homogeneous or zero Neumann boundary conditions. It will be shown in the application section that this is not a substantial restriction. Of coarse, it is important to use a suitable model which fulfills the assumed criteria as well as possible. For the sake of completeness a possible discretization of the element matrices arising from the boundary integrals over $\mathcal{A}_{N1}$ and $\mathcal{A}_{N2}$ are given in the appendix in Chapter B and Chapter C.

With the above considerations about the boundary integral terms the base equations used for the further finite element assembling can be written

$$\int_{\mathcal{V}}\left(\vec{\nabla}\times\vec{N}_i\right) \cdot\... ... - \jmath\omega\int_{\mathcal{V}}\mu\vec{N}_i\cdot\vec{\nabla}\psi \mathrm{d}V = 0$$ (5.33)

$$\int_{\mathcal{V}}\vec{\nabla}\lambda_i\cdot\left(\mu\vec{H}_1\ri... ...hcal{V}}\vec{\nabla}\lambda_i\cdot\left(\mu\vec{\nabla}\psi\right) \mathrm{d}V = 0.$$ (5.34)

$\vec{H}_1$ and $\psi$ are substituted by their approximations (5.27) and (5.28) in (5.33) and (5.34)

$$\quad \quad \sum_{j=1}^{m}c_j\int_{\mathcal{V}}\left(\vec{\nabla}... ...ega\sum_{j=1}^{m}c_j\int_{\mathcal{V}}\mu\vec{N}_i\cdot\vec{N}_j \mathrm{d}V -$$ (5.35)

$$\quad - \jmath\omega\sum_{j=m+1}^{n}c_j\int_{\mathcal{V}}\mu\vec{N}_i\cdot\vec{\nabla}\lambda_j \mathrm{d}V =$$

$$= - \sum_{j=n+1}^{M}c_j\int_{\mathcal{V}}\left(\vec{\nabla}\times... ...a\sum_{j=n+1}^{M}c_j\int_{\mathcal{V}}\mu\vec{N}_i\cdot\vec{N}_j \mathrm{d}V +$$

$$\quad + \jmath\omega\sum_{j=M+1}^{N}c_j\int_{\mathcal{V}}\mu\vec{N}_i\cdot\vec{\nabla}\lambda_j \mathrm{d}V , i\in[1;m]$$

$$\quad \quad \sum_{j=1}^{m}c_j\int_{\mathcal{V}}\vec{\nabla}\lambd... ...}\vec{\nabla}\lambda_i\cdot\left(\mu\vec{\nabla}\lambda_j\right) \mathrm{d}V =$$ (5.36)

$$= - \sum_{j=n+1}^{M}c_j\int_{\mathcal{V}}\vec{\nabla}\lambda_i\cd... ...\cdot\left(\mu\vec{\nabla}\lambda_j\right) \mathrm{d}V , i\in[m{+}1;n].$$

The linear equation system consisting of (5.35) and (5.36) is written more conveniently in matrix form

$$\begin{array}{cc} \left[ \begin{array}{ll} \left[A\right] & \... ...ray} \right] & \left\{c\right\} \end{array} = \left\{b\right\},$$ (5.37)

where (5.36) is multiplied by $-{\jmath}\omega$ to obtain a symmetric matrix. Thus, $\left[A\right]$ , $\left[B\right]$ , $\left[C\right]$ and $\left[b\right]$ are given by the expressions

$$A_{ij} = \int_{\mathcal{V}}\left(\vec{\nabla}\times\vec{N}_i\right) \cdo... ...t_{\mathcal{V}}\mu\vec{N}_i\cdot\vec{N}_j \mathrm{d}V, i\in[1;m], j\in[1;m]$$ (5.38)

$$B_{ij} = - \jmath\omega\int_{\mathcal{V}}\mu\vec{N}_i\cdot\vec{\nabla}\lambda_j \mathrm{d}V, i\in[1;m], j\in[m{+}1;n]$$ (5.39)

$$C_{ij} = \jmath\omega\int_{\mathcal{V}}\vec{\nabla}\lambda_i\cdot\left(\mu\vec{\nabla}\lambda_j\right) \mathrm{d}V, i\in[m{+}1;n], j\in[m{+}1;n]$$ (5.40)

$$b_i = - \sum_{j=n+1}^{M}c_j\int_{\mathcal{V}}\left(\vec{\nabla}\times... ...a\sum_{j=n+1}^{M}c_j\int_{\mathcal{V}}\mu\vec{N}_i\cdot\vec{N}_j \mathrm{d}V +$$ (5.41)

$$\quad + \jmath\omega\sum_{j=M+1}^{N}c_j\int_{\mathcal{V}}\mu\vec{N}_i\cdot\vec{\nabla}\lambda_j \mathrm{d}V , i\in[1;m]$$

$$b_i = \jmath\omega\sum_{j=n+1}^{M}c_j\int_{\mathcal{V}}\vec{\nabla}\lambda_i\cdot\left(\mu\vec{N}_j\right) \mathrm{d}V -$$ (5.42)

$$\quad - \jmath\omega\sum_{j=M+1}^{N}c_j\int_{\mathcal{V}}\vec{\na... ..._i\cdot\left(\mu\vec{\nabla}\lambda_j\right) \mathrm{d}V , i\in[m{+}1;n].$$

Obviously the right hand side vector $\{b\}$ can be calculated by (5.37) where $[A]$ , $[B]$ , and $[C]$ are expressed as in (5.38), (5.39), and (5.40) and the corresponding ranges for the global indexes $i$ and $j$ are taken from (5.41) and (5.42).