The Finite Element Method Applied to the Magnetostatic and Magnetodynamic Problems

Abstract

Modelling of realistic electromagnetic problems is presented by partial differential equations (FDEs) that link the magnetic and electric fields and their sources. Thus, the direct application of the analytic method to realistic electromagnetic problems is challenging, especially when modeling structures with complex geometry and/or magnetic parts. In order to overcome this drawback, there are a lot of numerical techniques available (e.g. the finite element method or the finite difference method) for the resolution of these PDEs. Amongst these methods, the finite element method has become the most common technique for magnetostatic and magnetodynamic problems.

1. Introduction

Mathematical modeling of realistic problems in the framework of electromagnetics leads to a set of partial derivates equations that have to be solved on a domain with complex geometry associated with boundary conditions and initial conditions. This complexity makes any analytical approach unpracticable. In the past (until 1960), people used experimentation (very expensive, sometimes destructive) or analogic simulation (lack of generality) to solve these problems. Since 1970, the growth of computer capabilities makes the numerical simulation a tool that is more and more used by the people interested in solving these complex problems. When using the computer, the continuous problem is represented with a finite number of degrees of freedom (d.o.f.). The continuous problem is then replaced by a discrete problem. There are a lot of numerical techniques available. We will see that the most common ones can be derived from the same general principle of weighted residuals.

A continuous formulation of a problem cannot generally be solved analytically and some numerical methods have to be used in order to obtain quantitative information about the solution. The unknown functions of a continuous problem belong to continuous function spaces which are usually of infinite dimensions, that is, those functions are usually described by an infinite number of parameters. The basis of any numerical method is to discretize such a problem in order to obtain a similar discrete problem, characterized by a finite number of unknowns which are called degrees of freedom. This discretization process consists of replacing the considered continuous function spaces by some discrete function spaces, whose dimensions are finite, and which are usually subspaces of them. Those spaces are also called approximation spaces and their elements are called approximation functions.

The function spaces are defined in a particular studied domain. If this one is discretized, that is, if it is defined as the union of geometric elements of simple shapes, and if the discrete function spaces are built in such a way that their functions are piecewise defined, then the approximation numerical method is called the finite element method (FEM). It is this kind of method we are interested in. We can thus see that the finite element method necessitates a double discretization: a discretization of some function spaces and a discretization of the studied geometric domain, which leads to a mesh.

Weak formulations are well adapted to the finite element method, which will appear in the following. Such formulations make use of several kinds of Green formulas.

2. Numerical technique

2.1 The Laplacian problem

The formalism used in the case of a Laplacian problem is sufficiently simple to be very understandable without lack of generality. The description of a Laplacian problem is presented now. Let us consider a bounded domain $\mathrm{\Omega }$ and its boundary $\mathrm{\Gamma }={\mathrm{\Gamma }}_{\mathrm{h}}$ ∪ ${\mathrm{\Gamma }}_{\mathrm{e}}$ (Figure 1).

Figure 1.

Studied domain $\mathrm{\Omega }$ and its boundary $\mathrm{\Gamma }$ = ${\mathrm{\Gamma }}_h\cup {\mathrm{\Gamma }}_e.$

The Laplace equation has to be solved in $\mathrm{\Omega }$ [1, 2, 3]:

$$\mathrm{\Delta }u\left(\mathbf{x}\right)=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}=0,$$E1

where u is the unknown field defined at each point x (x, y, z) of the studied domain. The associated boundary conditions are respectively Dirichlet and Neumann conditions, that is

$$u\left(\mathbf{x}\right)=\overline{u}\left(\mathbf{x}\right),\mathbf{x}\in {\mathrm{\Gamma }}_h,$$E2

$$v\left(\mathbf{x}\right)=\frac{\partial u\left(\mathbf{x}\right)}{\partial n}=\overline{v}\left(\mathbf{x}\right),\mathbf{x}\in {\mathrm{\Gamma }}_e.$$E3

This diffusion equation describes a wide range of physical phenomena. The next table shows some of these phenomena.

Problems	u	on ${\mathrm{\Gamma }}_h$	on ${\mathrm{\Gamma }}_e$
Thermostatics	T (temperature)	T fixed	thermal flux fixed
Electrostatics	V (voltage)	V fixed (electrode)	fixed flux of electrical displacement
Perfect fluids	$\psi$ (flow function); I (current function)
Magnetostatics	Br (reduced potential)	fixed magnetic flux density	fixed tangential magnetic field

A natural way to discretize the problem is to impose the error on the equation and on the boundary conditions weighted by a trial function w to be equal to zero, that is

$$\underset{\mathrm{\Omega }}{\int }\mathrm{\Delta }uwd\mathrm{\Omega }+\underset{{\mathrm{\Gamma }}_h}{\int }\left(\overline{u}-u\right)wd{\mathrm{\Gamma }}_e+\underset{{\mathrm{\Gamma }}_e}{\int }\left(\overline{v}-v\right)wd{\mathrm{\Gamma }}_h=0.$$E4

Equations (1–3) are then meanly solved, the sense of the mean being the principle of the numerical method. In fact, the numerical method used (F.D.M, F.E.M or B.E.M) are directly related to the chosen trial functions.

2.2 Green formulas

The following notations are used for integration of products of scalar or vector fields over a volume $\mathrm{\Omega }$ or on a surface $\mathrm{\Gamma }$, where L² and L² are the spaces of square-summable scalar and vector functions [2, 3]:

$$\left(\mathrm{u},\mathrm{v}\right)={\int }_{\mathrm{\Omega }}\mathrm{u}\left(\mathit{x}\right)\mathrm{v}\left(\mathit{x}\right)\mathrm{d}\mathit{x},\mathrm{u},\mathrm{v}\in {\mathrm{L}}^2\left(\mathrm{\Omega }\right)$$

$$\left(\mathit{u},\mathit{v}\right)={\int }_{\mathrm{\Omega }}\mathit{u}\left(\mathit{x}\right)·\mathit{v}\left(\mathit{x}\right)\mathrm{d}x,\mathit{u},\mathit{v}\in {\mathit{L}}^2\left(\mathrm{\Omega }\right),$$

$${⟨\mathrm{u},\mathrm{v}⟩}_{\mathrm{\Gamma }}={\int }_{\mathrm{\Gamma }}\mathrm{u}\left(\mathbf{x}\right)·\mathrm{v}\left(\mathbf{x}\right)\mathrm{ds},{⟨\mathbf{u},\mathbf{v}⟩}_{\mathrm{\Gamma }}={\int }_{\mathrm{\Gamma }}\mathrm{u}\left(\mathbf{x}\right)·\mathrm{v}\left(\mathbf{x}\right)\mathrm{ds},$$

A first relation of vectorial analysis

$$\mathbf{u}\bullet \text{gradv}+\mathrm{v}\bullet div\mathbf{u}=div\left(\mathrm{v}\mathbf{u}\right),$$

integrated in the domain $\mathrm{\Omega }$, after applying the divergence theorem, gives the Green formula said of kind grad-div in $\mathrm{\Omega }$, that is

$$\left(\mathit{u},gradv\right)+\left(div\mathit{u},v\right)=<v,\mathit{n}·\mathit{u}{>}_{\mathrm{\Gamma }},\forall \mathit{u}\in {\mathit{H}}^1\left(\mathrm{\Omega }\right),\forall v\in H^1\left(\mathrm{\Omega }\right)$$E5

where ${\mathit{H}}^1\left(\mathrm{\Omega }\right)\text{and}v\in H^1\left(\mathrm{\Omega }\right)$ are function spaces built for scalar and vector fields, respectively.

Another relation of vectorial analysis

$$\mathbf{u}.curlv-curl\mathit{u}·\mathit{v}=div\left(\mathit{v}\times \mathit{u}\right)$$E6

integrated in the domain $\mathrm{\Omega }$, after applying the divergence theorem, gives the Green formula said of kind curl-curl in $\mathrm{\Omega }$, that is

$$\left(\mathbf{u},curl\mathit{v}\right)-\left(curl\mathit{u},\mathit{v}\right)=<\mathit{u}\times \mathit{n},\mathit{v}{>}_{\mathrm{\Gamma }},.\forall \mathit{u}\mathbf{,}\mathit{v}\in {\mathit{H}}^1\left(\mathrm{\Omega }\right)$$E7

Note that the surface integral term appearing in this last formula can take the following similar forms:

$$⟨\mathit{u}\times \mathit{n},\mathit{v}{>}_{\mathrm{\Gamma }}={⟨\mathit{v}\times \mathit{u},\mathit{n}⟩}_{\mathrm{\Gamma }}=-⟨\mathit{v}\times \mathit{n},\mathit{u}{>}_{\mathrm{\Gamma }}$$

It is possible to define a generalized Green formula by

$$\left(\mathrm{Lu},\mathrm{v}\right)-\left(\mathrm{u},{\mathrm{L}}^{\ast }\mathrm{v}\right)={\int }_{\mathrm{\Gamma }}\mathrm{Q}\left(\mathrm{u},\mathrm{v}\right)\mathrm{ds},\forall \mathrm{u}\in dom\left(\mathrm{L}\right)\text{and}\forall \mathrm{v}\in dom\left({\mathrm{L}}^{\ast }\right),$$E8

where L and L* are differential operators of order n which act respectively on functions u and v defined in $\overline{\mathrm{\Omega }}$, with $\overline{\mathrm{\Omega }}$ = $\mathrm{\Omega }\cup \mathrm{\Gamma }$; Q is a bi-linear function of u and v. The operator L* is called the dual operator of L. It can easily be seen that formulas (6) and (7) are particular cases of (8).

2.3 Weak formulations

Consider a partial differential problem of the form [4].

$$\mathrm{L}\mathrm{u}=\mathrm{f}\text{in}\mathrm{\Omega },$$E9

$$\mathrm{B}\mathrm{u}=\mathrm{g}\mathrm{on}\mathrm{\Gamma }=\mathrm{\partial \Omega },$$E10

where L is a differential operator of order n, B is an operator which defines a boundary condition, f and g are functions respectively defined in Ω and on its boundary Γ, and u is an unknown function from a function space U and defined in $\overline{\mathrm{\Omega }}$, that is, u ∈ U($\overline{\mathrm{\Omega }}$). Note that f can eventually depend on u.

Problems (9 and 10) constitute what is called a classical formulation, or strong formulation. A function u ∈ U($\overline{\mathrm{\Omega }}$) which verifies this problem is called a classical solution, or strong solution. Particularly, as L is of order n, the function u has to be n–1 times continuously differentiable, that is, u ∈ Cn–1(Ω).

A weak formulation of problem (9) is defined as having the generalized form.

$$\left(\mathrm{u},{\mathrm{L}}^{\ast }\mathrm{v}\right)-\left(\mathrm{f},\mathrm{v}\right)+{\int }_{\mathrm{\Gamma }}{\mathrm{Q}}_{\mathrm{g}}\left(\mathrm{v}\right)\mathrm{ds}=0,\forall \mathrm{v}\in \mathrm{V}\left(\mathrm{\Omega }\right)$$E11

where L* is the dual operator of L, defined by the generalized Green formula (8), Qg is a linear form in v which depend on g, and the space V (Ω) is a space of test functions which has to be defined according to the operator L* and particularly according to the boundary condition (9 and 10). A function u which satisfies this equation for any test function v ∈ V (Ω) is called a weak solution.

The generalized Green formula (8) can be applied to formulation (11) in order to get L instead of L*, which usually consists of performing an integration by parts. It is then possible to find again, thanks to a judicious choice of test functions, the equations and relations of the classical formulation of the problem, that is, Eq. (9) and boundary condition (11).

It is often easy to check that a classical solution is also a weak solution. Nevertheless, it is not always straightforward that a weak solution is also a classical solution because it has to be regular enough in order to be defined at the classic sense.

One mathematical advantage of weak formulations is that they usually enable to prove the existence of a solution easier than classical formulations do. The solution has then to be proved to be regular enough to be also a classical solution. Another advantage of weak formulations is that they are well adapted to a discretization using finite elements and then to a numerical solution, which is not the case with classical formulations.

In some cases, it is possible to define a minimization problem similar to the weak formulation (11).

2.4 A weak formulation for the magnetodynamic problem

In order to illustrate the notion of weak formulation, consider the magnetodynamic problem, limited to the domain Ω, with boundary $\mathrm{\partial \Omega }$ = $\mathrm{\Gamma }$ = $\mathrm{\Gamma }$_h$\cup \mathrm{\Gamma }$_e (Figure 2), whose equations and material relations are written in Euclidean space ${\mathbb{R}}^3$[5, 6].

$$curl\mathbf{h}={\mathbf{j}}_s$$E12

$$curl\mathbf{e}=-\mathit{j\omega }\mathit{b}$$E13

$$div\mathit{b}=0$$E14

with behavior relations of materials.

Figure 2.

Studied domain Ω and its boundary.

$$\mathit{b}=\mu \mathit{h}$$E15

$$\mathit{j}=\sigma \mathit{e},$$E16

where $j$ is called the imaginary unit, $\mathit{b}$ is the magnetic induction (T), $\mathit{e}$ is the electric field (V/m), ${\mathit{j}}_s$ is the current density (A/m2), $\mathit{h}$ is the magnetic field (A/m), $\mu$ is the relative permeability and $\sigma$ is the electric conductivity (S/m). From the Eq. (13), the field $\mathit{b}$ can be obtained from a magnetic vector potential ${\mathit{a}}_i$ via the term:

$$\mathit{b}=curl\mathit{a}.$$E17

Combining (15 and 16) into (14), one has curl ($\mathit{e}+{\partial }_{\mathrm{t}}\mathit{a})=0$, that leads to the presentation of an electric scalar potential $\nu$ through $\mathit{e}=-{\partial }_{t}a-grad\upsilon$.

By starting from the Ampere’s law (12), the weak form of magnetic vector potential is written as [4, 6].

$${\left({\mathrm{\mu }}^{-1}\mathrm{b},curl{\mathrm{a}}^{\prime }\right)}_{\mathrm{\Omega }}-{\left(\mathrm{\sigma e},{\mathrm{a}}^{\prime }\right)}_{{\mathrm{\Omega }}_{\mathrm{c}}}+{⟨\mathrm{n}\times \mathrm{b},{\mathrm{a}}^{\prime }⟩}_{\mathrm{\Gamma }}={\left({\mathrm{j}}_{\mathrm{s}},{\mathrm{a}}^{\prime }\right)}_{{\mathrm{\Omega }}_{\mathrm{s}}}\forall {\mathrm{a}}^{\prime }\in {\mathrm{H}}_{\mathrm{e}}^0\left(curl,\mathrm{\Omega }\right)$$E18

Combining the magnetic vector potential $\mathit{a}$ and the electrical field $\mathit{e}$, one has

$${\left({\mu }^{-1}curl\mathit{a},curl{\mathit{a}}^{\prime }\right)}_{\mathrm{\Omega }}+{\left(\sigma {\partial }_{\mathrm{t}}\mathit{a},{\mathit{a}}^{\prime }\right)}_{{\mathrm{\Omega }}_c}+{\left(\sigma grad\upsilon ,curl{\mathit{a}}^{\prime }\right)}_{{\mathrm{\Omega }}_c}+{⟨\mathit{n}\times \mathit{h},{\mathit{a}}^{\prime }⟩}_{{\mathrm{\Gamma }}_h}={\left({\mathit{j}}_s,{\mathit{a}}^{\prime }\right)}_{{\mathrm{\Omega }}_s},\forall {\mathit{a}}_i^{\prime }\in F_e^0\left(curl,{\mathrm{\Omega }}_i\right),$$E19

where $H_e^0\left(curl,\mathrm{\Omega }\right)$is a function space defined on $\mathrm{\Omega }$ containing the basis functions for $\mathit{a}$ as well as for the test function ${\mathit{a}}^{\prime }$ (at the discrete level, this space is defined by edge FEs; notations (·, ·) and < ·, · > are respectively a volume integral in and a surface integral of the product of their vector field arguments.

2.5 A weak formulation for the magnetostatic problem

The magnetostatic problem is considered as a simplification of the magnetodynamic formulation where all time dependent phenomena are neglected. In a same way, by starting from the Ampere’s law (12), this initial form of the problem is its classical formulation.

Consider the Green formula of type grad-div in Ω (5) applied to the field b and to a scalar field ϕ’ to be defined, that is [6, 7, 8]

$$\left(\mathit{b},grad\mathrm{\varphi }\prime \right)+\left(div\mathit{b},\mathrm{\varphi }\prime \right)=⟨\mathit{n}\mathit{\bullet }\mathit{b},\mathrm{\varphi }\prime {>}_{\mathrm{\Gamma }},\forall \mathrm{\varphi }\prime \in \mathrm{\Phi }\left(\mathrm{\Omega }\right).$$E20

If the space Φ(Ω) is defined such as

$$\mathrm{\Phi }\left(\mathrm{\Omega }\right)=\left\{\mathrm{\varphi }\in \mathrm{H}1\left(\mathrm{\Omega }\right);\mathrm{\varphi }\mid {\mathrm{\Gamma }}_{\mathrm{h}}=0\right\},$$E21

then the last term of Eq. (20) is reduced to < n.b, ϕ’ > Γ_e and is equal to zero if condition (15) is taken into account. Moreover, the second term of this equation is equal to zero because of Eq. (16). Eq. (20) can then be reduced to.

$$\left(\mathit{b},grad\mathrm{\varphi }’\right)=0,\forall \mathrm{\varphi }’\in \mathrm{\Phi }\left(\mathrm{\Omega }\right).$$E22

This last form is called a weak formulation of the problem. It has been established starting from a Green formula but it can be considered now as an a priori posed form whose enclosed information can be deduced.

In fact, weak formulation (22) contains both Eq. (12) and boundarycondition (15). Indeed, by applying the Green formula of type grad-div to it, we get.

$$\left(div\mathit{b},\mathrm{\varphi }’\right)=⟨\mathit{n}\bullet \mathit{b},\mathrm{\varphi }\prime {>}_{\mathrm{\Gamma }},\forall \mathrm{\varphi }’\in \mathrm{\Phi }\left(\mathrm{\Omega }\right).$$E23

This equation is verified for any test function ϕ’ ∈ Φ(Ω) and thus, particularly, for any function ϕ’ whose value is equal to zero on Γ, that is, ϕ’ ∈ Φ₀ (Ω) because Φ₀(Ω) ⊂ Φ(Ω). Therefore, it comes that (div b, ϕ’) = 0 for any function ϕ’ of this kind and, consequently, that div b = 0 in Ω, that is, Eq. (12) is satisfied. Then, Eq. (23) is reduced to $⟨\mathit{n}\bullet \mathit{b},\mathrm{\varphi }\prime {>}_{\mathrm{\Gamma }}$= 0, and by considering now all the functions ϕ’ ∈ Φ(Ω) without any restriction, that is, which can vary freely on ${\mathrm{\Gamma }}_e$, it comes that ${\left.\mathit{n}\bullet \mathit{b}\right|}_{{\mathrm{\Gamma }}_e}$= 0, that is, that condition (13) is satisfied.

It is possible to obtain more information from the weak formulation, particularly as far as the transmission conditions on surfaces inside Ω are concerned. Consider for that two subdomains Ω₁ and Ω₂ of Ω separated by an interface Σ (Figure 3) [7].

Figure 3.

Interface Σ between two subdomains Ω₁ and Ω₂.

Let us apply the Green formula of type grad-div (5) to the fields b and ϕ’ successively in both subdomains Ω₁ and Ω₂. By taking into account that div b = 0 in Ω, and thus particularly in Ω₁ and Ω₂, then by summing the obtained relations, we get the relation [6, 7].

$$\left(\mathit{b},grad\mathrm{\varphi }’\right){\mathrm{\Omega }}_1\approx {\mathrm{\Omega }}_2=<\mathit{n}.\left({\mathit{b}}_{\mathit{1}}–{\mathit{b}}_{\mathit{2}}\right),\mathrm{\varphi }’>\mathrm{\Sigma }+<\mathit{n}.\mathit{b},\mathrm{\varphi }’>({\mathrm{\Omega }}_1\approx {\mathrm{\Omega }}_2),\forall \mathrm{\varphi }’\in \mathrm{\Phi }\left(\mathrm{\Omega }\right),$$E24

where b₁ and b₂ represent the field b on both sides of Σ in the respective domains Ω1 and Ω2. Considering the test functions ϕ’ whose support is Ω₁∪Ω₂ and which are equal to zero on _(Ω₁∪Ω₂), it remains from (24) the well known transmission condition n.(b₁–b₂)∣Σ = 0. Note that the first term of (24) vanishes thanks to Eq. (22) indeed, the domain of integration Ω₁ ∪ Ω₂ can be extended to Ω thanks to the chosen test functions.

The way to establish a weak formulation of Eq. (13) has been described and the richness of the information enclosed in such a formulation has been shown up. Using a similar procedure, a weak formulation associated with Eq. (12) can be established, but we will proceed differently in order to keep some classical equations.

If the field h is decomposed into a given source component h_s, such as curl h_s = j, and a reaction component h_r, then curl h_r = 0 and h_r is therefore of the form h_r = − grad ϕ (if Ω is simply connected). This consists of satisfying Eq. (15) classically. Taking into account the behavior law (15), we can write b = μ (hs − grad ϕ) and put this last expression in (24) to obtain.

$$\left(\mathrm{\mu }\left({\mathit{h}}_s-grad\mathrm{\varphi }\right),grad\mathrm{\varphi }’\right)=0,\forall \mathrm{\varphi }’\in \mathrm{\Phi }\left(\mathrm{\Omega }\right).$$E25

This formulation contains the whole problem (12 and 13). The potential ϕ is the unknown and all the other fields can be deduced from ϕ thanks to the equations which have been kept on a classical form. It appears that the potential ϕ belongs to the same space of the test functions or at least to a space Φ_r (Ω) which is parallel to it, that is, where the boundary condition relative to ϕ on Γh is not necessarily homogeneous, that is, ϕ∣Γ_h = constant. Note that this boundary condition on Γ_h is implicitly taken into account in the space Φ(Ω).

Weak formulation (25) can be considered as a system of an infinite number of equations with an infinite number of unknowns. It will be seen in the following how such a problem can be approximated to lead to a numerical solution. This approximation will constitute the phase of discretization.

A similar minimization problem

It is possible to define a minimization problem associated with (25). For that, let us define the functional [2, 3].

$$\mathrm{W}\left(\mathrm{\varphi }\right)=\left(\mathrm{\mu }\left(\mathbf{h}\mathrm{s}-grad\mathrm{\varphi }\right),{\mathit{h}}_s-grad\mathrm{\varphi }\right),$$E26

and let us pose the following minimization problem:

$$\text{find}\mathrm{\varphi }\in {\mathrm{\Phi }}_{\mathrm{r}}\left(\mathrm{\Omega }\right)\text{such}\mathrm{as}\mathrm{W}\left(\mathrm{\varphi }\right)\le \mathrm{W}\left(\mathrm{\varphi }’\right),\forall \mathrm{\varphi }’\in \mathrm{\Phi r}\left(\mathrm{\Omega }\right).$$E27

The physical materials are considered having linear magnetic behavior, but the following can be generalized easily for nonlinear materials.

By stationarizing functional (25) in relation to ϕ, it can be verified that (25) is obtained. It can also be verified that the solution ϕ of (25) minimizes this functional. Indeed, let us suppose that ϕ ∈ Φ_r(Ω) is solution of (25) and let us consider any function ψ ∈ Φ_r(Ω); then let us pose η = ψ − ϕ, which implies η ∈ Φ(Ω); we have

$$\mathrm{W}\left(\mathrm{\psi }\right)=\mathrm{W}\left(\mathrm{\varphi }+\mathrm{\eta }\right)=\left(\mathrm{\mu }\left({\mathit{h}}_s-grad\left(\mathrm{\varphi }+\mathrm{\eta }\right)\right),{\mathit{h}}_s-grad\left(\mathrm{\varphi }+\mathrm{\eta }\right)\right).$$

and thus.

$$\mathrm{W}\left(\mathrm{\psi }\right)=\mathrm{W}\left(\mathrm{\varphi }\right)+\left(\mathrm{\mu }grad\mathrm{\eta },grad\mathrm{\eta }\right)+\left(\mathrm{\mu }\left({\mathit{h}}_s-grad\mathrm{\varphi }\right),-grad\mathrm{\eta }\right).$$

As the second term of this sum is positive or equal to zero and the third term is equal to zero, because of (25), it comes that W(ψ) and W(ϕ).

Formulations (25) and (27) are then similar. Note that W(ϕ) is the magnetic coenergy and that the problem actually consists of minimizing this coenergy.

If the continuous function spaces are replaced by discrete spaces, and if the considered test functions are limited to these spaces, then the information inside a weak formulation will only be satisfied approximately, or weakly.

The basis of the discretization of weak formulations can be illustrated for the above magnetostatic problem, whose weak formulation is (25), that is

$$\left(\mathrm{\mu }\left({\mathit{h}}_s-grad\mathrm{\varphi }\right),grad\mathrm{\varphi }’\right)=0,\forall \mathrm{\varphi }’\in \mathrm{\Phi }\left(\mathrm{\Omega }\right),$$E28

with ϕ ∈ Φ(Ω). The space Φ(Ω) has to be replaced by a discrete space Φ_h(Ω) which is a subset of it, that is, Φ_h(Ω) ⊂ Φ(Ω). This space has a finite dimension, denoted N, and can then be defined by N linearly independent base functions. The principle is then to look for the function ϕ in Φ_h (Ω), which consists of determining N unknown parameters. This function will be only an approximation of the exact solution ϕ ∈ Φ (Ω). The more the functions of Φ_h (Ω) approximate well those of Φ (Ω), the higher the quality of the approximation is. Each test function ϕ’ will lead to an equation of the form (28) and, as the number of equations and unknowns has to be the same, N linearly independent test functions have to be chosen. This choice can be made on the base functions of Φ_h (Ω) and the method is called the Galerkine method. Such base functions are defined thanks to finite elements.

3. Finite elements

3.1 Definition of a finite element

A finite element is defined by the three element set (K, PK, ΣK) where [2, 6, 7]:

• K is a domain of space called a geometric element (usually of simple shape, that is, a tetrahedron, a hexahedron or a prism);

• PK is a function space of dimension nK, defined in K with base functions;

• ΣK is a set of nK degrees of freedom represented by nK linear functionals ϕi, 1 ≤ i ≤ nK, defined in space PK;

moreover, any function u ∈ PK must be defined uniquely by the degrees of freedom of ΣK, which defines the unisolvance of the finite element (K, PK, ΣK).

The role of a finite element is to interpolate a field in a function space of finite dimension. Several finite elements can be defined on the same geometric element and then, under certain conditions, can form mixed finite elements. Figure 4 shows the various spaces which occur in the definition of a finite element; the definition of the subspace of points κ ⊂ K is actually associated with the definition of the functionals.

Figure 4.

Spaces associated with a finite element (K, P_K, Σ_K) [6].

For the most commonly used finite elements, the degrees of freedom are associated with nodes of K and the functionals ϕi are reduced to functions of the coordinates in K; these elements are called nodal finite elements. Nevertheless, the above definition is more general thanks to the freedom let in the choice of the functionals. It will be shown that these can be, in addition to nodal values, integrals along segments, on surfaces or in volumes; the subspace of points κ ⊂ K (Figure 4) is then respectively a point, a segment, a surface or a volume.

3.2 Unisolvant finite element

The finite element (K, PK, ΣK) is unisolvant if [6].

$$\forall \mathrm{p}\in {\mathrm{P}}_{\mathrm{K}},\mathrm{\varphi i}\left(\mathrm{p}\right)=0;\forall \mathrm{\varphi i}\in \mathrm{\Sigma K}\Rightarrow \mathrm{p}\equiv 0.$$

In this case, for any function u regular enough, one can define a unique interpolation uK, called P_K-interpolant, such as.

$$\mathrm{\varphi i}\left(\mathrm{u}-{\mathsf{u}}_{\mathsf{K}}\right)=0,\forall \mathrm{\varphi i}\in {\mathrm{\Sigma }}_{\mathrm{K}};{\mathsf{u}}_{\mathsf{K}}\in {\mathrm{P}}_{\mathrm{K}}.$$E29

The set Σ_K is said P_K - unisolvant.

Proof:

Each function p ∈ P_K can be written as a linear combination of functions of a base of P_K, denoted {pi, 1 ≤ i ≤ n_K}, that is

$$\mathrm{p}=\sum _{\mathrm{i}=1}^{{\mathrm{n}}_{\mathrm{K}}}{\mathrm{a}}_{\mathrm{i}}{\mathrm{p}}_{\mathrm{i}},$$

where the pi, 1 ≤ i ≤ nK, are called base functions.

As the functionals ϕj, 1 ≤ j ≤ n_K, are linear, we have

$${\mathrm{\varphi }}_{\mathrm{j}}\left(\mathrm{p}\right)=\sum _{\mathrm{i}=1}^{{\mathrm{n}}_{\mathrm{K}}}{\mathrm{a}}_{\mathrm{i}}{\mathrm{\varphi }}_{\mathrm{j}}\left({\mathrm{p}}_{\mathrm{i}}\right),1\le \mathrm{j}\le {\mathrm{n}}_{\mathrm{K}}.$$

And, as ϕj(p) = 0, 1 ≤ j ≤ n_K, leads to p ≡ 0, the determinant of the matrix Φ (Φ_ji = ϕ_j(p_j), 1 ≤ i, j ≤ n_K) is not equal to zero; indeed the solution of the corresponding system must be identically equal to zero (i.e., ai = 0, 1 ≤ i ≤ n_K). Consequently, the system

$${\mathrm{\varphi }}_{\mathrm{j}}\left(\mathrm{u}\right)={\mathrm{\varphi }}_{\mathrm{j}}\left({\mathrm{u}}_{\mathrm{K}}\right)\iff {\mathrm{\varphi }}_{\mathrm{j}}\left(\mathrm{u}\right)=\sum _{\mathrm{i}=1}^{{\mathrm{n}}_{\mathrm{K}}}{\mathrm{a}}_{\mathrm{i}}{\mathrm{\varphi }}_{\mathrm{j}}\left({\mathrm{p}}_{\mathrm{i}}\right),1\le \mathrm{j}\le {\mathrm{n}}_{\mathrm{K}.}$$

has a unique solution (a_j, 1 ≤ i ≤ n_K).

3.3 Degrees of freedom

The interpolation of a function u, in the space P_K and in K, is given by the expression [4]

$${\mathrm{u}}_{\mathrm{K}}=\sum _{\mathrm{i}=1}^{{\mathrm{n}}_{\mathrm{K}}}{\mathrm{a}}_{\mathrm{i}}{\mathrm{p}}_{\mathrm{i}},{\mathrm{u}}_{\mathrm{K}}\in {\mathrm{P}}_{\mathrm{K}},$$

where the n_K coefficients a_j associated with the base functions p_j ∈ P_K can be determined thanks to relations (26), that is, thanks to the solution of the linear system.

$${\mathrm{\varphi }}_{\mathrm{j}}\left(\mathrm{u}\right)=\sum _{\mathrm{i}=1}^{{\mathrm{n}}_{\mathrm{K}}}{\mathrm{a}}_{\mathrm{i}}{\mathrm{\varphi }}_{\mathrm{j}}\left({\mathrm{p}}_{\mathrm{i}}\right),1\le \mathrm{j}\le {\mathrm{n}}_{\mathrm{K}},$$

provided that the function u is sufficiently regular for the ϕj(u), 1 ≤ j ≤ n_K, to exist.

This solution is simplified to the maximum if we define the functionals so that

$${\mathrm{\varphi }}_{\mathrm{j}}\left({\mathrm{p}}_{\mathrm{j}}\right)={\delta }_{\mathit{ij}},1\le \mathrm{i},\mathrm{j}\le {\mathrm{n}}_{\mathrm{K}}$$E30

where δ_i,j is the Kronecker symbol, that is

$${\mathrm{\delta }}_{\mathrm{ij}}=\left\{\begin{array}{c}1s_{i}i=j\\ 0s_{i}i\ne j\end{array}\right.$$

The matrix of the system is then the unit matrix and the solution is

$${\mathrm{a}}_{\mathrm{j}}={\mathrm{\varphi }}_{\mathrm{j}}\left(\mathrm{u}\right),1\le \mathrm{j}\le {\mathrm{n}}_{\mathrm{K}}$$

In this case, the interpolation u_K ∈ P_K is expressed by

$${\mathrm{u}}_{\mathrm{K}}=\sum _{\mathrm{j}=1}^{{\mathrm{n}}_{\mathrm{K}}}{\mathrm{\varphi }}_{\mathrm{j}}\left(\mathrm{u}\right){\mathrm{p}}_{\mathrm{j}},$$E31

where the coefficients ϕj(u) = ϕj(u_K), 1 ≤ j ≤ n_K, are called degrees of freedom.

3.4 Finite element spaces

A finite element space X_h can be built on a set of geometric elements and associated finite elements. Its definition depends on the mesh M_h of the domain Ω as well as the knowledge of the finite element (K, P_K, Σ_K) associated with each domain K ∈ M_h [6]

Given a function u defined in Ω, regular enough, its interpolant u_h ∈ X_h is uniquely defined such as [6]:

• The restriction ${\left.{\mathrm{u}}_{\mathrm{h}}\right|}_{\mathrm{K}}$, that is, the form of u_h in the geometric element K, belongs to the space P_K;

• The restriction Finite element spaces

A finite element space X_h can be built on a set of geometric elements and associated finite elements. Its definition depends on the mesh M_h of the domain Ω as well as the knowledge of the finite element (K, P_K, Σ_K) associated with each domain K∈ M_h

Given a function u defined in Ω, regular enough, its interpolant u_h ∈ X_h is uniquely defined such as [6]:

• The restriction ${\left.{\mathrm{u}}_{\mathrm{h}}\right|}_{\mathrm{K}}$, that is, the form of uh in the geometric element K, belongs to the space P_K;

• The restriction ${\left.{\mathrm{u}}_{\mathrm{h}}\right|}_{\mathrm{K}}$ is entirely determined by the knowledge of the set of values Σ_K(u) of the degrees of freedom of the function u - this is a consequence of the unisolvance;

Some continuous conditions have to be ensured across the interfaces between geometric elements, which is the property of conformity.

A mesh M_h of the studied domain Ω is defined as a collection of geometric elements which have in common either a facet, or an edge, or a node, or nothing (Figure 5). The elements cannot overlap each other.

Figure 5.

Mesh of a part of a two-dimensional domain Ω.

The finite element space X_h has a finite dimension, denoted Dh. It can be characterized by a set of degrees of freedom Σ_h linked up to the sets Σ_K,∀K∈ M_h, that is

$${\mathrm{\Sigma }}_{\mathrm{h}}=\left\{{\mathrm{\varphi }}_{\mathrm{h},\mathrm{j}},1\le \mathrm{j}\le {\mathrm{D}}_{\mathrm{h}}\right\}.$$

It is also possible to define the base functions p_h,j, 1 ≤ i ≤ Dh, of the space Xh from the base functions of the spaces P_K, ∀ K∈ Mh. Those have to verify the relations

$${\mathrm{\varphi }}_{\mathrm{j}}\left({\mathrm{p}}_{\mathrm{j}}\right)={\delta }_{\mathit{ij}},1\le \mathrm{i},\mathrm{j}\le {\mathrm{n}}_{\mathrm{K}}$$E32

similar to relations (32). They are actually piecewise defined and their supports are as “small” as possible, that is, are constituted by a limited number of geometric elements.

Then, with any function u regular enough so that the degrees of freedom ϕ_h,j(u), 1 ≤ j ≤ D_h, are well defined, it can be associated a function u_h, called X_h-interpolant, defined by

$${\mathrm{u}}_{\mathrm{h}}=\sum _{\mathrm{j}=1}^{{\mathrm{D}}_{\mathrm{h}}}{\mathrm{\varphi }}_{\mathrm{h},\mathrm{j}}\left(\mathrm{u}\right){\mathrm{p}}_{\mathrm{h},\mathrm{j}}$$E33

4. Construction of a sequence of finite element spaces

4.1 Geometric elements

A mesh of a domain is considered which is built with a collection of geometric elements which can be tetrahedra (4 nodes), hexahedra (8 nodes) and prisms (6 nodes) (Figure 6) [4, 6, 9].

Figure 6.

Collection of different geometric elements [6].

These elements are called volumes and their vertices represent nodes. The sets of nodes, edges, facets and volumes of this mesh are denoted by N, E, F and V, respectively. Their sizes are #N, #E, #F and #V.

The i-th node of the mesh is denoted by ni or {i}. The edges and facets can be defined with ordered sets of nodes. An edge is denoted by eij or {i, j}, a triangular facet by fijk or {i, j, k}, and a quadrangular facet by fijkl or {i, j, k, l}. These geometric entities are shown in Figure 7.

Figure 7.

Geometric entities: Node, edge and facets (i, j, k, l ∈ N) [6].

4.2 Base functions of spaces Sⁱ

Consider the function pi(x) of coordinates of point x and relative to node ni, which is equal to 1 at this node, varies continuously in geometric elements having this node in common, and becomes equal to 0 in other elements without discontinuity (Figure 6). This function is nothing else than the base function, relative to node ni, of the function space of nodal finite elements built on the considered geometric elements. The function subspaces associated with each of the finite elements have respective dimensions 4, 8 or 6, for tetrahedra, hexahedra and prisms [4, 6, 9].

With node n_i = {i}, is associated the function

$${\mathrm{s}}_{{\mathrm{n}}_{\mathrm{i}}}\left(\mathrm{x}\right)={\mathrm{p}}_{\mathrm{i}}\left(\mathrm{x}\right).$$E34

The finite dimensional space generated by all sn_i’s is denoted by S⁰.

With edge e_ij = {i, j}, is associated the vector field

$${\mathrm{s}}_{{\mathrm{e}}_{\mathrm{ij}}}={\mathrm{p}}_{\mathrm{j}}grad\sum _{\mathrm{r}\in {\mathrm{N}}_{\mathrm{F},\mathrm{j}\overline{\mathrm{i}}}}{\mathrm{p}}_{\mathrm{r}}-{\mathrm{p}}_{\mathrm{i}}grad\sum _{\mathrm{r}\in {\mathrm{N}}_{\mathrm{F},\mathrm{i}\overline{\mathrm{j}}}}{\mathrm{p}}_{\mathrm{r}},$$E35

where ${\mathrm{N}}_{\mathrm{F},\mathrm{m}\overline{\mathrm{n}}}$ is the set of nodes which belong to the facet of the geometrical element including evaluation point x, and including node m but not node n; such a facet is uniquely defined for three-edge-per-node elements. Its determination is shown in Figure 8, where either a triangular or a quadrangular facet is involved, and where shown edges belong to the geometric element including point x. Directions of dotted edges can be modified in order to schematize either a tetrahedron, a hexahedron or a prism. The defined set of nodes comes into view as being either {{m}, {o}, {p}} or {{m}, {o}, {p}, {q}}, respectively. The set ${\mathrm{N}}_{\mathrm{F},\mathrm{m}\overline{\mathrm{n}}}$ depends on point x, thus on elements. Particularly, it is empty (no node) in elements which have not edge {m, n} in common. Consequently, field s_eij is zero in all the elements non adjacent to edge e_ij.

Figure 8.

Determination of the facet associated with ${\mathrm{N}}_{\mathrm{F},\mathrm{m}\overline{\mathrm{n}}}$ [6].

The vector field space generated by all s_e is denoted by S¹.

With facet f = f_ijk = {i, j, k} = {q₁, q₂, q₃} or f = f_ijkl = {i, j, k, l} = {q₁, q₂, q₃, q₄}, is associated the vector field

$${\mathrm{s}}_{\mathrm{f}}={\mathrm{a}}_{\mathrm{f}}\sum _{\mathrm{c}=1}^{\#{\mathrm{N}}_{\mathrm{f}}}{\mathrm{p}}_{{\mathrm{q}}_{\mathrm{c}}}grad\left(\sum _{\mathrm{r}\in {\mathrm{N}}_{\mathrm{F},{\mathrm{q}}_{\mathrm{c}}}\overline{{\mathrm{q}}_{\mathrm{c}}}+1}{\mathrm{p}}_{\mathrm{r}}\right)\times grad\left(\sum _{\mathrm{r}\in {\mathrm{N}}_{\mathrm{F},{\mathrm{q}}_{\mathrm{c}}}\overline{{\mathrm{q}}_{\mathrm{c}}}-1}{\mathrm{p}}_{\mathrm{r}}\right)$$E36

where #Nf is the number of nodes of facet f, a_f = 2 if #N_f = 3, a_f = 1 if #N_f = 4, and the list of q_i’s is made circular by setting q₀ ≡ q#N_f and q#N_f + 1 ≡ q₁. Field s_f is zero in all the elements non adjacent to facet f.

Vector fields s_fijk(l)'s generate the space S².

With volume v, is associated the function sv, equal to 1/vol(v) on v and 0 elsewhere. The space S3 is generated by these functions.

Some developments give the following results: sn_i is equal to 1 at node ni, and to 0 at other nodes; the circulation of s_eij is equal to 1 along edge eij, and to 0 along other edges; the flux of s_fijk(l) is equal to 1 across facet s_fijk(l), and to 0 across other facets; and the volume integration of sv is equal to 1 over volume v, and to 0 over other volumes; that is

$${\mathrm{s}}_{\mathrm{i}}\left({\mathrm{x}}_{\mathrm{j}}\right)={\mathrm{\delta }}_{\mathrm{ij}},\forall \mathrm{i},\mathrm{j}\in \mathrm{N}$$E37

$${\int }_{\mathrm{j}}{\mathrm{s}}_{\mathrm{i}}·\mathrm{dl}={\mathrm{\delta }}_{\mathrm{ij}},\forall \mathrm{i},\mathrm{j}\in \mathrm{E}$$E38

$${\int }_{\mathrm{j}}{\mathrm{s}}_{\mathrm{i}}·\mathbf{n}\mathrm{ds}={\mathrm{\delta }}_{\mathrm{ij}},\forall \mathrm{i},\mathrm{j}\in \mathrm{F}$$E39

$${\int }_{\mathrm{j}}{\mathrm{s}}_{\mathrm{i}}\mathrm{dv}={\mathrm{\delta }}_{\mathrm{ij}},\forall \mathrm{i},\mathrm{j}\in \mathrm{V}$$E40

where δ_ij = 1 if i = j and δ_ij = 0 if i ≠ j.

These properties show up various kinds of functionals and involve that functions s_n, s_e, s_f, s_v form bases for the spaces they generate. They are then called nodal, edge, facet and volume base functions. The associated finite elements are called nodal, edge, facet and volume finite elements.

4.3 Geometric interpretation of edge and facet functions

A geometric interpretation of edge and facet functions may be helpful to verify some of their properties. The vector field [4, 6, 9]

$$grad{\mathrm{P}}_{\mathrm{F},\mathrm{m}\overline{\mathrm{n}}}=grad\sum _{\mathrm{r}\in {\mathrm{N}}_{\mathrm{F},\mathrm{m}\overline{\mathrm{n}}}}{\mathrm{p}}_{\mathrm{r}},$$E41

involved in both expressions (31) and (32), should be analyzed at first. The continuous scalar field,

$${\mathrm{P}}_{\mathrm{F},\mathrm{m}\overline{\mathrm{n}}}=\sum _{\mathrm{r}\in {\mathrm{N}}_{\mathrm{F},\mathrm{m}\overline{\mathrm{n}}}}{\mathrm{p}}_{\mathrm{r}},$$E42

has the characteristic of being equal to 1 at every point on the facet associated with ${\mathrm{N}}_{\mathrm{F},\mathrm{m}\overline{\mathrm{n}}}$. This is a property of the nodal base functions. Therefore, vector field (42) is orthogonal to this facet at every point on it (Figure 9).

Figure 9.

Geometric interpretation of the edge function s_e(35) [6].

The vector field which is the product of pm and (35),

$${\mathrm{p}}_{\mathrm{m}}grad\sum _{\mathrm{r}\in {\mathrm{N}}_{\mathrm{F},\mathrm{m}\overline{\mathrm{n}}}}{\mathrm{p}}_{\mathrm{r}},$$E43

is considered now. This field is said to be associated with edge {m, n}. As far as the function pm is concerned, it is equal to 0 on all the edges of the geometric element including point x, except those which are incident to node {m}. Therefore, the circulation of (43) is equal to 0 along all the edges except emn; field (43) is either simply equal to zero on them, or orthogonal to them (Figure 9). The combination of two fields of form (43) associated with edges {j, i} and {i, j}, as in (35), leads to a vector field which has the same properties as (43) (Figure 9), and has consequently the announced properties of s_eij. The fact that its circulation along edge eij is equal to 1 needs some calculation to be proved.

The vector field

$${\mathrm{p}}_{{\mathrm{q}}_{\mathrm{c}}}grad{\mathrm{P}}_{\mathrm{F},{\mathrm{q}}_{\mathrm{c}}}{\overline{\mathrm{q}}}_{\mathrm{c}+1}\times grad{\mathrm{P}}_{\mathrm{F},{\mathrm{q}}_{\mathrm{c}}{\overline{\mathrm{q}}}_{\mathrm{c}-1}}$$E44

which appears in expression (35) of s_f, is considered now. Both gradients in (44) are shown in Figure 10. Each one is orthogonal to its associated facet and, therefore, their cross product (i.e., a × b in Figure 10) is parallel to both these facets.

Figure 10.

Vector field a × b involved in s_f(33) [6].

The flux of this cross product, and in consequence the one of (44), is then equal to 0 across these facets. The term pq_c in (44) enables the flux of (44) to be equal to zero across all other facets except facet f. The summation in (44) keeps the same property. The flux of s_f across facet f is then the only one to differ from zero (Figure 11).

Figure 11.

Geometric interpretation of the facet function s_f(36) [6].

4.4 Degrees of freedom

The expression of a field in the base of a space Sⁱ –S⁰ or S³ for a scalar field, S¹ or S² for a vector field– gives scalar coefficients, called degrees of freedom. Fields ϕ ∈ S⁰, h ∈ S¹, j ∈ S² and ρ ∈ S³ can be expressed as [4, 6, 9]

$$\mathrm{\Phi }=\sum _{\mathrm{n}\in \mathrm{N}}{\mathrm{\varphi }}_{\mathrm{n}}{\mathrm{s}}_{\mathrm{n}},\mathrm{\varphi }\in {\mathrm{S}}^0,{\mathrm{\varphi }}_{\mathrm{n}}=\mathrm{\varphi }\left(x_{\mathrm{n}}\right),\mathrm{n}\in \mathrm{N},$$E45

$$\mathbf{h}=\sum _{\mathrm{e}\in \mathrm{E}}{\mathrm{h}}_{\mathrm{e}}{\mathrm{s}}_{\mathrm{e}},\mathbf{h}\in {\mathrm{S}}^1,{\mathrm{h}}_{\mathrm{e}}={\int }_{\mathrm{e}}\mathbf{h}·\mathrm{d}\mathbf{l},\mathrm{e}\in \mathrm{E}$$E46

$$\mathbf{j}=\sum _{f\in F}{j}_f{\mathrm{s}}_{\mathbf{f}},\mathbf{j}\in S^2,{\mathrm{j}}_{\mathrm{f}}={\int }_{\mathrm{f}}\mathbf{j}·\mathbf{n}\mathrm{ds},\mathrm{f}\in \mathrm{F}$$E47

$$\mathrm{\sigma }=\sum _{\mathrm{v}\in \mathrm{V}}{\mathrm{\sigma }}_{\mathrm{v}}{\mathrm{s}}_{\mathrm{v}},\mathrm{\rho }\in {\mathrm{S}}^3,{\mathrm{\sigma }}_{\mathrm{v}}={\int }_{\mathrm{v}}\mathrm{\sigma dv},\mathrm{v}\in \mathrm{V}$$E48

The degrees of freedom ϕ_n, h_e, j_f and ρ_v are thus, respectively, values at nodes, circulations along edges, fluxes across facets or volume integrals, of the associated fields. This is a consequence of the base functions. The associated linear functionals, as mentioned in the definition of finite elements, are thus respectively pointwise evaluations, line, surface and volume integrals.

4.5 Continuity of base functions across facets

It can be proved that the function sn is continuous across facets. The same holds true for the tangential component of s_e and for the normal component of s_f. As for function sv, it is discontinuous. This property, called conformity, allows to take exactly into account interface conditions for fields used in the modeling of physical problems. For example, in electromagnetic problems, vector fields of S¹ can represent vector fields like magnetic field h or electric field e whose tangential components are continuous across interfaces between materials, and those of S² can represent fields like induction field b or current density field j whose normal components are continuous across interfaces between these materials.

4.6 Spaces Sⁱ form a sequence

The notion of incidence is first defined [4, 6, 9]:

The incidence of node n in edge e, denoted by i (n, e), is equal to 1 if n is the extremity of e, −1 if n is the origin of e, and 0 if n does not belong to e.

Next, the incidence of edge e in facet f is denoted by i(e, f). If e belongs to f, and if the ordered set of nodes of e appears as a direct subset in the circular set of nodes of f, then it is equal to 1. It is equal to −1 in the case of an inverse subset. If e does not belong to f, it is equal to 0.

Finally, the incidence of facet f in volume v is denoted by i(f, v). If f belongs to v, and if the normal to f, whose direction is given by the ordered set of nodes of f (right-hand rule), is outer to v, then it is equal to 1. It is equal to −1 in the case of an inner normal. If f does not belong to v, it is equal to 0.

Thanks to this notion, the following equalities can be proved,

$$\sum _{\mathrm{e}\in \mathrm{E}}\mathrm{i}\left(\mathrm{n},\mathrm{e}\right){\mathrm{s}}_{\mathrm{e}}=grad{\mathrm{s}}_{\mathrm{n}}$$E49

$$\sum _{\mathrm{f}\in \mathrm{F}}\mathrm{i}\left(\mathrm{e},\mathrm{f}\right){\mathrm{s}}_{\mathrm{f}}=curl{\mathrm{s}}_{\mathrm{e}},$$E50

$$\sum _{\mathrm{v}\in \mathrm{V}}\mathrm{i}\left(\mathrm{f},\mathrm{v}\right){\mathrm{s}}_{\mathrm{v}}=div{\mathrm{s}}_{\mathrm{f}}.$$E51

The following inclusions are then verified,

$$grad\left({\mathrm{S}}^0\right)\subset {\mathrm{S}}^1,curl\left({\mathrm{S}}^1\right)\subset {\mathrm{S}}^2,div\left({\mathrm{S}}^2\right)\subset {\mathrm{S}}^3.$$E52

Therefore, the spaces Si, i = 0 to 3, form a sequence, that can be schematized by the diagram in Figure 12.

Figure 12.

The sequence of spaces Sⁱ.

These spaces can then constitute approximation spaces for some continuous spaces Fi, i = 0 to 3, which contain scalar and vector fields associated with electromagnetic fields. The associated finite elements can then be called mixed elements.

All the established properties of base functions are valid for any collection of considered geometric elements, that is, for any mixing of tetrahedra, hexahedra and prisms.

5. Practical information about finite elements

5.1 Isoparametric elements

An isoparametric element is a finite element whose nodal base functions, which enable the interpolation of scalar fields, are also used to parametrize the associated geometric element. The base functions are usually piecewise defined, in each of the geometric elements which cover the studied domain, and some continuity conditions have to be satisfied at the interfaces between elements. Then, there will be no discontinuity of the interpolated scalar fields, nor of the coordinates after transformation from the reference elements towards the real ones. Such base functions are said to be conformal.

Consider a nodal finite element (K, PK, ΣK). If NK is the set of nodes of K, whose coordinates are xi, i ∈ N, and if the pi(u), i ∈ NK, are its base functions expressed in the coordinates u of the reference element Kr associated with K, then the parametrization of K (i.e., x = x(u)) is given by [6]

$$\mathbf{x}=\sum _{\mathrm{i}\in {\mathrm{N}}_{\mathrm{K}}}{\mathbf{x}}_{\mathrm{i}}{\mathrm{p}}_{\mathrm{i}}\left(\mathbf{u}\right)$$E53

where x∈K, u∈Kr; this element is isoparametric.

5.2 Reference elements

We define here the reference elements which are associated with the considered geometric elements, that is, with tetrahedra, hexahedra and prisms. Nodal, edge, facet and volume finite elements are defined in these geometric elements.

5.2.1 Reference tetrahedron of type I

The reference tetrahedron of type I is an element with 4 nodes whose coordinates are given in Figure 13. The associated geometric entities, as well as their notation, are shown in Figure 13. The nodal and edge base functions of this element are given in Tables 1 and 2. Table 3 shows the notation of facets. The incidence matrices are given by (53), (54) and (55) (Figure 14).

Figure 13.

Reference tetrahedron of type I [6].

Node; i∈N	Nodal base function; pi (u, v, w) = si (u, v, w)
1	1 − u − v − w
2	u
3	v
4	w

Table 1.

Nodal base functions of the tetrahedron of type I.

Edge	e = {i, j}		se(u), u = (u, v, w)
e∈E	i∈N	j∈N	se,u	se,v	se,w
1	1	2	1 − v –w	u	u
2	1	3	v	1 − u –w	v
3	1	4	w	w	1 − u –v
4	2	3	– v	u	0
5	2	4	– w	0	u
6	3	4	0	– w	v

Table 2.

Notation of the edges of the tetrahedron of type I and associated edge base functions (s_e).

Facet	f = {i, j, k}
f∈F	i∈N	j∈N	k∈N
1	1	2	4
2	1	3	2
3	1	4	3
4	2	3	4

Table 3.

Notation of the facets of the tetrahedron of type I.

Figure 14.

Geometric entities defined on a tetrahedron of type I [6].

Edge-node incidence matrix

E54

Facet-edge incidence matrix

E55

Volume-facet incidence matrix

E56

5.2.2 Reference hexahedron of type I

Figure 15.

Reference hexahedron of type I [6].

Figure 16.

Geometric entities defined on a hexahedron of type I [4].

Node; i∈N	Nodal base function; pi (u, v, w) = si (u, v, w)
1	(1 − u) (1 − v) (1 − w) / 8
2	(1 + u) (1 − v) (1 − w) / 8
3	(1 + u) (1 + v) (1 − w) / 8
4	(1 − u) (1 + v) (1 − w) / 8
5	(1 − u) (1 − v) (1 + w) / 8
6	(1 + u) (1 − v) (1 + w) / 8
7	(1 + u) (1 + v) (1 + w) / 8
8	(1 − u) (1 + v) (1 + w) / 8

Table 4.

Nodal base functions of the hexahedron of type I.

Edge	e = {i, j}			se(u), u = (u, v, w)
e∈E	i∈N	j∈N		se,u	se,v	se,w
1	1	2		(1 − v) (1 − w) / 8	0	0
2	1	4		0	(1 − u) (1 − w) / 8	0
3	1	5		0	0	(1 − u) (1 − v) / 8
4	2	3		0	(1 + u) (1 − w) / 8	0
5	2	6		0	0	(1 + u) (1 − v) / 8
6	3	4		–(1 + v) (1 − w) / 8	0	0
7	3	7		0	0	(1 + u) (1 + v) / 8
8	4	8		0	0	(1 − u) (1 + v) / 8
9	5	6		(1 − v) (1 + w) / 8	0	0
10	5	8		0	(1 − u) (1 + w) / 8	0
11	6	7		0	(1 + u) (1 + w) / 8	0
12	7	8		–(1 + v) (1 + w) / 8	0	0

Table 5.

Notation of the edges of the hexahedron of type I and associated edge base functions (s_e).

Facet	f = {i, j, k, l}
f∈F	i∈N	j∈N	k∈N	l∈N
1	1	2	6	5
2	1	4	3	2
3	1	5	8	4
4	2	3	7	6
5	3	4	8	7
6	5	6	7	8

Table 6.

Notation of the facets of the hexahedron of type I.

The reference hexahedron of type I is an element with 8 nodes whose coordinates are given in Figure 15. The associated geometric entities, as well as their notation, are shown in Figure 16. The nodal and edge base functions of this element are given in Tables 4 and 5. Table 6 shows the notation of facets. The incidence matrices are given by (56), (57) and (58).

Edge-node incidence matrix

E57

Facet-edge incidence matrix

E58

Volume-facet incidence matrix

E59

5.2.3 Reference prism of type I

The reference prism of type I is an element with 6 nodes whose coordinates are given in Figure 17. The associated geometric entities, as well as their notation, are shown in Figure 18. The nodal and edge base functions of this element are given in Tables 7 and 8. Table 9 shows the notation of facets. The incidence matrices are given by (59), (60) and (61).

Figure 17.

Reference prism of type I [6].

Figure 18.

Geometric entities defined on a prism of type I [6].

Node; i∈N	Nodal base function; pi (u, v, w) = si (u, v, w)
1	(1 − u − v) (1 − w) / 2
2	u (1 − w) / 2
3	v (1 − w) / 2
4	(1 − u − v) (1 + w) / 2
5	u (1 + w) / 2
6	v (1 + w) / 2

Table 7.

Nodal base functions of the prism of type I.

Edge	e = {i, j}			se(u), u = (u, v, w)
e∈E	i∈N	j∈N		se,u	se,v	se,w
1	1	2		(1 − v) (1 − w) / 2	u (1 − w) / 2	0
2	1	3		v (1 − w) / 2	(1 − u) (1 − w) / 2	0
3	1	4		0	0	(1 − u − v) / 2
4	2	3		– v (1 − w) / 2	u (1 − w) / 2	0
5	2	5		0	0	u / 2
6	3	6		0	0	v / 2
7	4	5		(1 − v) (1 + w) / 2	u (1 + w) / 2	0
8	4	6		v (1 + w) / 2	(1 − u) (1 + w) / 2	0
9	5	6		– v (1 + w) / 2	u (1 + w) / 2	0

Table 8.

Notation of the edges of the prism of type I and associated edge base functions (s_e).

Facette	f = {i, j, k (, l)}
f∈F	i∈N	j∈N	k∈N	l∈N
1	1	2	5	4
2	1	3	2	—
3	1	4	6	3
4	2	3	6	5
5	4	5	6	—

Table 9.

Notation of the facets of the prism of type I.

Edge-node incidence matrix

E60

Facet-edge incidence matrix

E61

Volume-facet incidence matrix

E62

6. Applications

The practical test problem is a 3-D model based on the benchmark problem 19 of the TEAM workshop including a stranded inductor (coil) and an aluminum plate (Figure 19) [10].

Figure 19.

Modeling of TEAM problem 7: Coil and conducting plate [10].

The coil is excited by a sinusoidal current which generates the distribution of time varying magnetic fields around the coil (Figure 20). The relative permeability and electric conductivity of the plate are ${\mu }_{r,\mathit{\text{plate}}}=1$,${\sigma }_{r,\mathit{\text{plate}}}=35.26\mathrm{MS}/\mathrm{m}$, respectively. The source of the magnetic field is a sinusoidal current with the maximum ampere turn being 2742AT. The problem is tested with two cases of frequencies of the 50 Hz and 200 Hz.

Figure 20.

The 3-D mesh model with edge elements of the coil and conducting plate, and the limited boundary [4] (left), and distribution of magnetic flux density generated by the excited sinusoidal current in the coil, with ${\mu }_{r,\mathit{\text{plate}}}=1,{\sigma }_{r,\mathit{\text{plate}}}=35.26\frac{\mathit{MS}}{m}\mathit{\text{and}}$ f = 50 Hz.

The 3-D dimensional mesh with edge elements is depicted in Figure 21 (left). The distribution of magnetic flux density generated by the excited electric current in the coil is pointed out in Figure 21 (right). The computed results on the of the z-component of the magnetic flux density along the lines A1-B1 and A2-B2 (Figure 19) is checked to be close to the measured results for different frequencies of exciting currents (already proposed by authors in [10]) are shown in Figure 21. The mean errors between calculated and measured methods [10] on the magnetic flux density are lower than 10%.

Figure 21.

The comparison of the calculated results with the measured results on magnetic flux densities at y = 72 mm, with ${\mu }_{r,\mathit{\text{plate}}}=1,{\sigma }_{r,\mathit{\text{plate}}}=35.26\frac{\mathit{MS}}{m}\mathit{\text{and}}$ different frequencies [4].

The y component of the varying of the eddy current losses with different frequencies (50 Hz and 200 Hz) along the lines A3-B3 and A4-B4 (Figure 19) is shown in Figure 22. The computed results are also compared with the measured results as well [7]. The obtained results from the theory modeling are quite similar as what measured from the measurements. The maximum error near the end of the conductor plate on the eddy currents between two methods are below 20% for both cases (50 Hz and 200 Hz).

Figure 22.

The comparison of the calculated results with the measured results at z = 19 mm, with ${\mu }_{r,\mathit{\text{plate}}}=1,{\sigma }_{r,\mathit{\text{plate}}}=35.26\frac{\mathit{MS}}{m}\mathit{\text{and}}$ different frequencies [4].

7. Conclusions

In the 3D computation of the magnetic flux density and eddy current, thanks to the set of Maxwell’s equations, it has been successfully developed for two weak formulations, where the discretization of the fields is performed by Whitney edge elements [2, 3, 8]: magnetostatic formulation and magnetodynamic formulation. The developments of the method is validated on the actual problem (TEAM problem 7) [10]. The numerical error between simlated and measured results on the magnetic flux densities and eddy current is lower than 10%. This is also proved that there is a very good validation between two methods. The results have been achieved by a detailed study of the magnetodynamic formulation.