[PDF] - Sliding Interfaces for Eddy Current Simulations Raffael Casagrande PDF Document

SLIDE 1

Sliding Interfaces for Eddy Current Simulations

Raffael Casagrande

Master Thesis Supervisor:

Prof. Dr. Ralf Hiptmair

Zürich, April 2013

SLIDE 2

2.1. Coulomb gauged formulation . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2. Temporal gauged formulation . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3. H-formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4. The eddy current problem in a moving, solid body . . . . . . . . . . . . . 9

3. Discontinuous Galerkin formulation 11

3.1. Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.1.1. Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.1.2. Broken spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.1.3. Discrete Friedrichs inequality . . . . . . . . . . . . . . . . . . . . . 17 3.1.4. Trace inequalities and polynomial approximation . . . . . . . . . . 17 3.2. DG formulation of Eddy Current problem . . . . . . . . . . . . . . . . . . 18 3.2.1. Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

4. Aspects of the Implementation 32

4.1. Reduction to 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

SLIDE 3

CONTENTS 4.2. Details of the imlementation. . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.2.1. The mesh data structure . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3. Variational formulation for moving systems . . . . . . . . . . . . . . . . . 40

5. Results 45

5.1. Convergence to analytical solution . . . . . . . . . . . . . . . . . . . . . . 45 5.2. Translational motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.3. Rotational motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

6. Conclusion and Outlook 54

A. Maxwell’s equations in a moving frame of reference

55

B. Documentation of the code

60

C. Acknowledgements

61 Bibliography 62

ii

SLIDE 4

1. Introduction

We present a discontinuous Galerkin (DG) approach for solving the transient eddy current problem in the presence of moving, rigid bodies which share common interfaces. Such problems are important in the modelling of electromechanical devices. In particular the induced currents can have positive or negative effects and are often essential for the functioning of electric engines. Consider for example the setting depicted in Fig. 1.1. Here multiple bodies are moving with respect to each other and sliding interfaces

exist. We attach to each body a separate coordinate frame and mesh so we can solve

the equations in Lagrangian variables. Figure 1.1.: An example configuration with 3 solid bodies. Separate coordinate systems are attached to each one. In order to deal with the non-matching grids we use symmetric-interior penalty tech- niques at the interfaces. We use conforming finite element spaces in the interior of each

body. Together with the appropriate transformation formulas we then link the different

descriptions into one consistent framework. This approach doesn’t need any remeshing and allows us to model the problem from a Lagrangian perspective. In particular no convective terms appear in the formulation and the method is therefore expected to be stable for bodies moving at high speed. In time, we use an implicit Euler discretization which results in a sparse, symmetric- positive-definite linear system that has to be solved in every iteration. The proposed discretization is analyzed in a general 3 dimensional DG framework but only 2 dimensional model problems have been implemented and tested. In space, the problem is discretized by means of first order Lagrangian elements as well as first order Nedelec

SLIDE 5

1. Introduction

elements of the first kind. Another commonly used approach to solve the same physical problem is based on mortar- methods (see for instance [14],[15], [4],[8]). These methods have typically symmetric- positive definite matrices as well and need no remeshing. However the implementation is harder because additional degrees of freedom are introduced along the interface. This is especially true for translational settings where the interface itself is changing over

time. There exist many other approaches which rely on conforming meshes, but they

almost always need a partial or full remeshing in every time step which increases the computational cost heavily (cf. [14, Introduction]). More recently Ferrer [7],[1] has applied the DG theory to the Navier-Stokes equations for treating sliding interfaces. The material is presented in five chapters: We start off in chapter 2 by giving a brief introduction to the eddy current model and its possible formulations. The chapter ends with a discussion of the Galilean transformation properties of the previously presented eddy current formulations. Chapter 3 starts by introducing mathematical tools from the DG world and continues with a presentation of the DG variational formulation. The latter is then analyzed and a proof of convergence is given. Following this rather theoretical material we present some implementation details in chapter 4. Finally, chapter 5 presents and discusses all the numerical results. The thesis ends with a short summary and some concluding remarks (chapter 6). 2

SLIDE 6

2. Eddy Current Problem

In this chapter we will derive three different formulations of the transient eddy current problem by neglecting the displacement current in Maxwell’s equations. By requiring a unique solution we will derive appropriate initial and boundary conditions for the given

problems. In the last part of this chapter we will show that the three formulations are

invariant under rotation and translation if the quantities are transformed correctly. It should be mentioned that there are many other ways of treating the eddy current

problem. In this short presentation, we focus our attention on the Coulomb gauged

potential formulation (unknowns: A, ϕ), the temporal gauged potential formulation (unknown: A) and the H formulation (unknown: the magnetic field H). We refer the reader to Albanese and Rubinacci [2] for a survey of different formulations.

Maxwell’s equations

The eddy current problem is a simplified version of Maxwell’s equations, Magnetic Gauss’s law: div B = 0 (2.1a) Faraday’s law of induction: curl E + ∂B ∂t = 0 (2.1b) Electric Gauss’s law: div D = ρ (2.1c) Ampère’s law: curl H − ∂D ∂t = jf + ji. (2.1d) Here B denotes the magnetic induction ([ Vs/m2]), E the electric field ([ V/m]), D the electric induction ([ N/Vm]), ρ the free charge density ([ C/m3]), H the magnetic field ([ A/m]), jf the free current density ([ A/m2]) and ji the impressed current density ([ A/m2]). The latter is a prescribed current density which can be used to model coils in the computational domain. We consider only linear, isotropic materials with immediate response in this thesis and therefore the constitutive relations D = εE B = µH jf = σE

SLIDE 7

2. Eddy Current Problem
apply. ε, µ and σ are all material dependent properties and are termed permittivity,

permeability and conductivity, respectively. They are simple scalars which are only a function of the material. In the following we are going to derive the eddy current problem. For this we need to make the central Assumption 2.1. The displacement current, ∂D

∂t in Ampères law (2.1d) is zero. This

is justified if the electric field is varying only slowly (Quasistatic simplification). This assumption is typically justified for very high conductivities respectively if the electric field E is varying only slowly.

2.1. Coulomb gauged formulation

Remark. In this chapter we assume that A ∈ [C∞(Ω)]3 and ϕ ∈ C∞(Ω) to simplify the

derivations and manipulations. This assumption is relaxed later when we consider the variational formulation in chapter 3. In addition we assume that the domain Ω is simply connected. By virtue of the Helmholtz decomposition (cf. [18]) we can express B and E through the vector potential A∗ and the scalar potential ϕ∗: B = curl A∗, (2.2a) E = − grad ϕ∗ − ∂A∗ ∂t . (2.2b) This way the first two Maxwell equations (2.1a),(2.1b) are fulfilled automatically and

nly Ampères law remains. Expressed in the variables A∗ and ϕ∗ it reads as

σ∂A∗ ∂t + curl 1 µ curl A∗ + σ grad ϕ∗ = ji. (2.3) It is well known that this equation doesn’t possess a unique solution. Indeed, we can use for example the gauge transformation A′∗ = A∗ + grad ψ, ϕ′∗ = ϕ∗ − ∂ψ ∂t , to transform A and ϕ without affecting the observable quantities B and E (ψ(x, t) is an arbitrary function). In practice the Coulomb gauge condition, div A∗ = 0 (2.4) is often introduced to fix this ambiguity. Together with appropriate boundary conditions this fixes the value of A and ϕ and leads to a well posed problem. We will refer to the system of equations represented by Eq. 2.3 together with Eq. 2.4 as the Coulomb gauged eddy current problem. 4

SLIDE 8

2.2. Temporal gauged formulation

Remark. If σ = 0 the Eqs. (2.3) and (2.4) become the elliptic boundary value problem

curl 1 µ curl A∗ = ji, div A∗ = 0. This shows that the eddy current model is principally not capable of modelling electromagnetic waves.

2.2. Temporal gauged formulation

Another possibility to guarantee that Eq. (2.3) has a unique solution consists in using the temporal gauge transformation A := A∗ +

t

grad ϕ∗(t)dt, ϕ := 0. (2.5) The electromagnetic fields can then be calculated through the relations B = curl A and E = −∂A ∂t . The gauge transformation 2.5 effectively corresponds to a gauge condition in which ϕ ≡ 0. If we apply the temporal gauge transform to the Coulomb gauged eddy current problem we arrive at the temporally gauged eddy current problem, curl 1 µ curl A + σ∂A ∂t = ji (2.6)

Boundary and Initial Conditions

So far we have omitted the boundary and initial conditions for the eddy current problem tactically. In this section we will derive such conditions by requiring Eq. (2.6) to have a unique solution. Notice that if A1 and A2 are both solutions of Eq. (2.6) then ˜ A = A1 − A2 must be a solution of the homogeneous equation, curl 1 µ curl ˜ A + σ∂ ˜ A ∂t = 0. Multiplying the equation by ˜ A, integrating over Ω and time and using Stokes theorem (cf. [20, Thm. 3.6]) we get

Ω

T

σ∂ ˜ A ∂t · ˜ A +

T

Ω

curl ˜ A · 1 µ curl ˜ A =

T

∂Ω

1 µ

curl ˜

A × n

· ˜

A. 5

SLIDE 9

2. Eddy Current Problem

The integral over time of the first term can be evaluated exactly if σ is independent of time:

Ω

1 2σ

˜

A(T)

2 − 1

2σ

˜

A(0)

2 + T

Ω

curl ˜ A · 1 µ curl ˜ A =

T

∂Ω

1 µ

curl ˜

A × n

· ˜

A. (2.7) From the above equation it is clear that the cases σ = 0 and σ > 0 must be handled

separately. For this we split the domain Ω into two distinct subdomains Ω = Ω1 ∪ Ω2

such that σ|Ω1 > 0 and σ|Ω2 ≡ 0. In the following we will therefore discuss the solution to (2.6) for each subdomain separately and consider the coupling of the two solutions later. In Ω1 we prescribe the following initial and boundary conditions: A(x, y, z, t = 0) = A(0)

n Ω1,

(2.8a) curl A × n = g1

n ∂Ω1,

(2.8b) This causes the second term on the left-hand side and the boundary term on the right- hand side of Eq. (2.7) to drop out:

Ω1

1 2σ

˜

A(T)

2 + T

Ω1

curl ˜ A · 1 µ curl ˜ A = 0 Because there are only quadratic terms left, ˜ A(T) ≡ 0 and curl ˜ A ≡ 0 in Ω1, i.e. A1 = A2. Therefore Eq. (2.6) together with the boundary/initial conditions (2.8) has a unique solution. In Ω2 things are more involved; Here we only require the boundary condition, curl A × n = g2

n ∂Ω2,

(2.9) so that Eq. (2.7) reduces to

T

Ω2 curl ˜

A · 1

µ curl ˜

A = 0. Once again the terms on the left hand side are quadratic, so that curl ˜ A

!

≡ 0 in Ω2. This implies that ˜ A is irrotational and therefore there exists a scalar function ψ such that ˜ A = grad ψ. To restore uniqueness one usually requires additionally the Coulomb gauge condition div A = 0

n Ω2 with b.c.

(2.10) A · n = h

n ∂Ω2.

(2.11) This fixes ψ to zero or a constant function so that ˜ A = 0 and therefore A1 = A2. It remains to connect the subproblems on Ω1 and Ω2 by suitable transmission conditions. For this we define Γ12 := Ω1 ∩ Ω2. Let p be a planar plane with boundary ∂p and define 6

SLIDE 10

2.2. Temporal gauged formulation l = Γ12 ∩ p to be the intersection of Γ12 and p. Now integrating Eq. (2.6) over p, assuming that no surface currents exist and using Stokes theorem we get: 0 =

p

(·) −

p∩Ω1

(·) −

p∩Ω2

(·) = −

l

1 µ|Ω1 curl A|Ω1 · dτ1 −

l

1 µ|Ω2 curl A|Ω2 · dτ2 = −

l
1

µ|Ω1 curl A|Ω1 − 1 µ|Ω2 curl A|Ω2

· dτ1

where τ1 = −τ2 denotes the line elements on either side of Γ12. Because p is arbitrary this results in the jump condition

1

µ curl A

× n12 = 0.

(2.12) where n12 points from Ω1 to Ω2. I.e. the tangential component of 1

µ curl A is continuous

across Γ12 (cf. Lemma 3.8 in chapter 3 for a more elegant and general proof of this jump/interface condition). From this and the boundary conditions (2.8b) and (2.9) we see that g1(x, y, z, t)|Γ12 = g2(x, y, z, t)|Γ12. We are now in the position to state the eddy current problem in temporal gauge with the appropriate boundary and initial conditions: curl 1 µ curl A + σ∂A ∂t = ji

n Ω

(2.13a) div A = 0

n Ω2

(2.13b) A(x, y, z, t = 0) = A(0)

n Ω

(2.13c) 1 µ curl A × n = g

n ∂Ω

(2.13d) A · n = h

n ∂Ω2

(2.13e)

1

µ curl A

× n12 = 0
n Γ12 = Ω1 ∩ Ω2

(2.13f) Remarks.

By taking the divergence of Eq. (2.13a) we see that div ji must fulfil the com-

patibility conditions div ji = 0 in Ω2. In addition h must fulfil the Neumann compatibility condition

∂Ω2 h ds = 0.
For proving the uniqueness of the IBVP (2.13) we assumed that σ is independent
f time.

This assumption is probably not needed to guarantee uniqueness but simplifies the proof.

This is the first possible interpretation of the eddy current equation. In contrast

to the interpretation presented next, this formulation works for σ = 0 as well as for σ > 0. 7

SLIDE 11

2. Eddy Current Problem

2.3. H-formulation

One can give a simpler but more restrictive meaning to the eddy current equation (2.13a) by dividing it on both sides by σ, applying the curl operator and making Assumption 2.2. We assume σ > 0 in Ω. This implies that Ω2 = ∅ and Ω1 = Ω. curl 1 σ curl 1 µ curl A + ∂ curl A ∂t = curl 1 µji, ⇔ curl 1 σ curl H + ∂µH ∂t = curl 1 σji

n Ω

(2.14) Of course we have to give a new meaning to the gauge condition (2.13b) and rephrase it in terms of H. By taking a glimpse at the magnetic Gauss law we easily state the new gauge condition: div µH = 0

n Ω.

(2.15) By using a similar reasoning as before for the temporally gauged formulation it is possible to show that the following initial-boundary value problem (IBVP) has a unique solution: curl 1 σ curl H + ∂µH ∂t = curl 1 σji

n Ω

(2.16a) div µH = 0

n Ω

(2.16b) H(x, y, z, t = 0) = H(0)

n Ω

(2.16c) H × n = g

n ∂Ω

(2.16d)

Remark. The gauge condition (2.16b) in this formulation is necessarily required for the

whole computational domain Ω by virtue of the magnetic gauss law. However by taking the divergence of Eq. (2.16a) we see that ∂ div µH

∂t

= 0. So if the gauge condition is fulfilled for time t = 0 it is automatically fulfilled for all t > 0. Usually one still includes the gauge condition in practice for the sake of numerical stability. It remains to show that the solution to the IBVP (2.16) fulfils Maxwell’s equations under assumptions 2.1 and 2.2. First we notice that Ampères law (2.1d) reduces to curl H = jf + ji = σE + ji. Solving this relation for E and inserting it into Faraday’s law of induction (2.1b) we get the eddy current equation (2.16a). I.e. Faraday’s law is fulfilled and jf = curl H − ji E = jf σ .

Remark. If we compare the continuity equation for electric charges, div (jf + ji) = ∂ρ

∂t ,

with the currents computed by the eddy current model we see ∂ρ

∂t = div curl H = 0. I.e.

the eddy current model preserves space charges. 8

SLIDE 12

2.4. The eddy current problem in a moving, solid body

Later on we are interested in settings where one solid body is moving inside the laboratory frame. It is not straightforward to calculate the electromagnetic fields in such a

setting. In particular Lorentz forces appear as source terms in the eddy current equa-
tions. For example the Coulomb gauged eddy current equation (2.13a) in the laboratory

frame of reference reads as (cf. appendix A, Eq. (A.7)) σ∂A∗ ∂t + curl 1 µ curl A∗ + σ grad ϕ∗ = ji + σV × curl A∗. (2.17) In appendix A we show that this equation is invariant under translation and rotation and takes exactly the same form in the moving frame of reference, σ∂ ˜ A∗ ∂t + ˜ curl 1 µ ˜ curl ˜ A∗ + σ ˜ grad ˜ ϕ∗ = ˜ ji + σ ˜ V × ˜ curl ˜ A∗, (2.18) where we used a tilde to denote field quantities which are defined in the moving frame

f reference. The field quantities of the laboratory frame are related to the fields in the

body’s coordinate frame by the relations T˜ E = E + V × B T ˜ B = B T˜ ji = ji T ˜ H = H T˜ jf = jf T ˜ A∗ = A∗ ˜ ϕ = ϕ − V · A∗ (2.19) Here T(t) stands for the rotation matrix of the moving frame (the principal coordinates in the moving system are the columns of T) and V stands for the velocity of of the moving frame (both are defined in the laboratory frame). Usually we will pick the moving frame of reference such that it coincides the body’s principal coordinate axes. In that case ˜ V ≡ 0, i.e. the Lorentz forces disappear in the moving frame of reference. In other words the Lorentz force term is due to a change in coordinates. The H-formulation is of course also invariant under rotation and translation and the quantities transform according to the transformation rules (2.19). It remains to derive a transformation formula for the temporal gauged formulation of the eddy current problem. The idea is to start from the Coulomb gauged formulation and use the temporal gauge transformation (2.5) to express A in the moving frame of reference. In other words we are imposing the gauge condition ˜ ϕ ≡ 0 in the moving frame. We can then use the transformation rules (2.19) to expand the expression in terms of the laboratory 9

SLIDE 13

2. Eddy Current Problem

coordinates: ˜ A := ˜ A∗ +

t

˜ grad ˜ ϕ(t)dt = TT A∗ +

t

TT grad ϕ −

t

TT grad (V · A∗) = TT A −

t

TT grad (V · A∗) (2.20) Where we have used the transformation rule for the del-operator, ˜ ∇ = TT ˜ ∇. We sum- marize the transformation rule in the Lemma 2.3 (Transformation rule for A). The temporally gauged eddy current problem with σ > 0, curl 1 µ curl A + σ∂A ∂t = ji + σV × curl A (2.21a) A(x, y, z, t = 0) = A(0) (2.21b) curl A × n = g (2.21c) is invariant under translation and rotation if ji transforms according to the rules (2.19) and if A and g transform as T ˜ A = A − T

t

TT grad (V · A). (2.22) T˜ g = g (2.23)

Proof. We can interpret A as a solution of (2.17) with A∗ = A, ϕ∗ = 0. Owing to

relation (2.20) transformation rule (2.22) follows immediately and it is clear that ˜ A is a solution to (2.21). The transformation rule for g follows from the relation T ˜ curl( TT g) = curl g (see Appendix A, Eq. (A.4)).

Remark. If we were to use the simple transformation rule T ˜

A = A instead of Eq. (2.22) the system would actually change gauge when switching from the laboratory to the moving frame of reference. In the rest frame the system is gauged such that ϕ ≡ 0 but in the moving frame ˜ ϕ = −V · ˜

A. I.e. an additional term would appear if the eddy

current equation was expressed in the principal coordinates of the moving frame. To

vercome this problem we applied the temporal gauge transformation to ˜

A and ˜ ϕ. We have thus provided transformation formulas for all quantities appearing in the Coulomb gauged, the temporal gauged and the H-formulation of the eddy current equa-

tion. In chapter 4 we will use these relations to derive the appropriate transmission

conditions at sliding interfaces. 10

SLIDE 14

3. Discontinuous Galerkin formulation

In the previous chapter we have seen that the solution of the eddy current problem is not necessarily continuous across interfaces. For instance in the Eddy current problem

nly the tangential components of the magnetic field H are continuous across interfaces.

The normal component of the field can jump. Thus edge/Nédelec elements (see [3], [13,

Sect. 5.5]) are often used in the context of finite element methods. They are tangentially

continuous by construction and are thus conforming in the space H(curl, Ω) (cf. Sect. 3.1.2). In the framework of discontinuous Galerkin methods one allows the numerical solution to be discontinuous across element faces. In particular we will allow the numerical solution of the eddy current problem to contain tangential discontinuities. Consequently this solution will not lie in the space H(curl, Ω) any more and thus the method is non-conforming. Because the numerical solution isn’t required to be continuous, additional conditions are required to assert uniqueness and convergence. There are various ways to achieve this; We will focus on the symmetric-weighted-interior-penalty (SWIP) method which essentially punishes tangential jumps of the numerical approximation. Using certain assumptions on the analytical solution and on the mesh we will be able to prove the convergence of the approximative solution to the analytical solution in a well defined sense (Sec. 3.2.1). The relaxed constraint on the continuity of the tangential components has the big advan- tage that non-matching meshes can easily be included in the discussion. This is because the resulting formulation is local to each element and the elements are interconnected

nly by their fluxes. As we will see later, the discontinuous Galerkin approach is indeed

very well suited for non-conforming meshes resulting from sliding interface problems.

3.1. Preliminaries

The discontinuous Galerkin approach relies on an extensive machinery of mathematical

tools. We will present many of them without proof in this section and refer the reader

to the excellent book of Di Pietro and Ern [6] for further information.

SLIDE 15

3. Discontinuous Galerkin formulation

3.1.1. Notation Later on we will need the following spaces: Lp(Ω) :=

f : Ω → R : f is Lebesque measurable,
Ω

|f(x)|p dµ(x) < ∞

,

W m,p(Ω) := {f ∈ Lp(Ω)| ∂αf ∈ Lp(Ω) for all |α|ℓ1 ≤ m} , Hm(Ω) := W m,2(Ω), H(curl, Ω) : =

f ∈
L2(Ω)

3

curl fL2(Ω) < ∞
where ∂αf and curl refer to the weak/distributional derivative of f and α is a multi-
index. The spaces are equipped with the norms and seminorms

fLp(Ω) :=

Ω

|f|p

1/p

fW m,p(Ω) :=

 

|α|ℓ1≤m

∂αfp

Lp(Ω)

 

1/p

|f|W m,p(Ω) :=

 

|α|ℓ1=m

∂αfp

Lp(Ω)

 

1/p

fHm(Ω) := fW m,2(Ω) |f|Hm(Ω) = |f|W m,2(Ω) fH(curl,Ω) :=

f2

L2(Ω) + |f|2 H(curl,Ω)

1/2

|f|H(curl,Ω) := curl fL2(Ω) . We used the convention that the L2 norm and the Hm norms of vector valued function f =

d

i=1 fi2. We also define the L2 inner product by (f, g) :=

Ω fg, and equiv-

alently for vector fields (f, g) :=

Ω f · g. For time dependent problems we consider the

space Cl(V ) := Cl([0, tF ]; V )

f V -valued functions which are l-times continuously differentiable in the interval [0, tF ].

I.e. if ψ ∈ Cl(V ), then ψ : (0, tF ) ∋ t → ψ(t) ≡ ψ(·, t) ∈ V.

Meshes

We will assume that there exists a mesh of the domain Ω and define it’s faces as (cf. Di Pietro and Ern [6, Definition 1.16]) Definition 3.1 (Mesh faces). Let Th be a mesh of the domain Ω. We say that a (closed) subset F in Ω is a mesh face if F has positive (d−1)-dimensional Hausdorff measure (in dimension 1, this means that F is nonempty) and if either one of the following conditions is satisfied:

1. There are distinct mesh elements T1 and T2 such that F = ∂T1 ∩ ∂T2; in such a

case, F is called an interface 12

SLIDE 16

3.1. Preliminaries

2. there is T ∈ Th such that F = ∂T ∩∂Ω; in such a case, F is called a boundary face.

Based on this we introduce the following sets of faces: Fi

h := set of interfaces in Th,

Fb

h := set of all boundary faces in Th,

Fh := Fi

h ∪ Fb h,

FT := {F ∈ Fh|F ⊂ ∂T} . For a given element/face T ∈ Th or F ∈ Fh we denote by hT resp. hF the diameter of the T resp. F. Further we define h := max

T∈Th ht.

Later on we will state polynomial approximation properties and trace inequalities which rely on a sequence of meshes TH := (Th)h∈H where the sequence H = (hn)n∈N is strictly decreasing and satisfies hn > 0 ∀n and hn → 0 as n → ∞. Further we assume that the mesh sequence TH is admissible (see [6, Definition 1.57]) and weakly quasi-uniform: Definition 3.2 (Weakly quasi-uniform mesh sequence (cf. [12])). A mesh sequence H is weakly quasi-uniform if there is a constant µ∗ with 0 < µ∗ < 1 such that hh−µ∗

min → 0 as h → 0

where hmin = min

T∈Th hT

One can obtain an admissible and weakly quasi-uniform sequence by starting out with a coarse, tetrahedral nonconforming mesh and then subdividing every tetrahedron uni- formly in 8 smaller tetrahedrons. In 2D the same procedure works with triangles. On each mesh element we can further define the polynomial space Pk

d(T),

Pk

d(T) :=

  p : Rd ⊃ T ∋ x → R

p(x) =
|α|ℓ1≤k

γαxα

  

where we used the multi-index α and T ∈ Th is a mesh element. Following the idea of Di Pietro and Ern [6] we define the broken polynomial space by Pk

d(Th) :=

v ∈ L2(Ω)
∀T ∈ Th, v|t ∈ Pk

d(T)

.

Later on in this chapter we will look for discrete solutions Ah to the eddy current problem (2.13) in broken polynomial space

Pk

3(Th)

3. This means that Ah can be discontinuous

across faces in tangential and normal direction. This motivates the introduction of a 13

SLIDE 17

3. Discontinuous Galerkin formulation

tangential jump resp. weighted average of a function Ah ∈

Pk

3(Th)

3 across a face

F ∈ Fh: for F ∈ Fi

h, F = ∂Ti ∩ ∂Tj :

for F ∈ Fb

h, F = ∂Ti ∩ ∂Ω :

[Ah]T = nF ×

Ah|Ti − Ah|Tj
[Ah]T = nF × Ah|Ti

(jump) for F ∈ Fi

h, F = ∂Ti ∩ ∂Tj :

for F ∈ Fb

h, F = ∂Ti ∩ ∂Ω :

{Ah}ω = ω1 Ah|Ti + ω2 Ah|Tj {Ah}ω = Ah|Ti . (average) where nF always points from Ti to Tj and ωi ∈ [0, 1] are weighting factors which fulfil the constraint ω1 + ω2 = 1. Notice that the tangential jump as well as the weighted average both act on vectors and produce vectors. In a similar fashion we define the usual jumps and averages for scalar valued functions ϕ ∈ Pk

3(Th),

for F ∈ Fi

h, F = ∂Ti ∩ ∂Tj :

for F ∈ Fb

h, F = ∂Ti ∩ ∂Ω :

[ϕ] = ϕ|Ti − ϕ|Tj [ϕ] = ϕ|Ti (jump) for F ∈ Fi

h, F = ∂Ti ∩ ∂Tj :

for F ∈ Fb

h, F = ∂Ti ∩ ∂Ω :

{ϕ}ω = ω1 ϕ|Ti + ω2 ϕ|Tj {ϕ}ω = ϕ|Ti . (average)

Remark. In the above definitions the corresponding traces of Ah, ϕ always exist because

Ah resp. ϕ lies in the broken polynomial space. Later on we will use the same definition

f the tangential jump and weighted average also for functions belonging to more general

function spaces. In order for this to work the traces of the function must exist. E.g. it doesn’t make sense to define a the tangential jump of a L2 function over an interface because there is no suitable definition for the trace of the function. We will fix this ambiguity by placing additional assumptions on the analytical solution and defer the details to Sec. 3.2.1.

Remark. For the scalar valued jumps it is necessary to fix the orientation of nF uniquely

for every edge. 3.1.2. Broken spaces In this section we extend the concepts of Sobolev spaces and the H(curl, Ω) space by allowing discontinuities across element faces (assuming there exists a suitable mesh). Following Di Pietro and Ern [6, Section 1.2.5] we define the broken Sobolev spaces as W m,p(Th) := {v ∈ Lp(Ω) | ∀T ∈ Th, v|T ∈ W m,p(T)} , Hm(Th) :=

v ∈ L2(Ω)
∀T ∈ Th, v|T ∈ Hm(T)
.

Analogous we define the broken version of the H(curl, Ω) space by H(curl, Th) :=

v ∈
L2(Ω)

3

∀T ∈ Th, v|T ∈ H(curl, T)
.

14

SLIDE 18

3.1. Preliminaries In order to simplify notation we introduce the broken gradient and broken curl operator. Definition 3.3 (Broken gradient, cf. [6]). The broken gradient grad h : W 1,p(Th) → [Lp(Ω)]d is defined such that, for all v ∈ W 1,p(Th), ∀T ∈ Th, ( grad hv)|T := grad (v|T ) . Definition 3.4 (Broken curl). The broken curl operator curlh : H(curl, Th) → [L2(Ω)]3 is defined such that, for all v ∈ H(curl, Th), ∀T ∈ Th, (curlh v)|T := curl (v|T ) . Intuitively it is clear that the Sobolev spaces W m,p(Ω) and Hm(Ω) are subspaces of the more general broken Sobolev spaces W m,p(Th) and Hm(Th). The same holds for the H(curl, Ω) spaces. The following statements are more precise and allow us to characterize the involved spaces. Lemma 3.5 (Broken gradient on usual Sobolev spaces). Let m ≥ 0 and let 1 ≤ p ≤ ∞. There holds W m,p(Ω) ⊂ W m,p(Th). Moreover, for all v ∈ W 1,p(Ω), grad hv = grad v in [Lp(Ω)]d. Lemma 3.6 (Characterization of W 1,p(Ω)). Let 1 ≤ p ≤ ∞. A function v ∈ W 1,p(Th) belongs to W 1,p(Ω) if and only if [v] := 0 ∀F ∈ Fi

h.

Lemma 3.7 (Broken curl operator on usual H(curl, Ω) space). There holds H(curl, Ω) ⊂ H(curl, Th). Moreover for all v ∈ H(curl, Ω), curlh v = curl v in [L2(Ω)]3.

Proof. Let v ∈ H(curl, Ω) and Φ ∈ [C∞

0 (T)]3 with T ∈ Th. Using the definition of the

broken curl operator and the notion of a weak derivative we can write:

T

curlh v · Φ =

T

curl (v|T ) · Φ =

T

v · curl Φ Let now EΦ be the extension of Φ by zero to Ω.

T

curlh v · Φ =

Ω

v · curl EΦ =

Ω

curl v · EΦ =

T

(curl v)|T · Φ I.e. curlh (v|T ) = (curl v)|T since Φ is arbitrary. Because T ∈ Th is arbitrary as well, curlh v = curl v. Lemma 3.8 (Characterization of H(curl, Ω)). A function v ∈ H(curl, Th)∩[W 1,1(Th)]3 belongs to H(curl, Ω) if and only if [v]T = 0 ∀F ∈ Fi

h

15

SLIDE 19

3. Discontinuous Galerkin formulation
Proof. Let v ∈ H(curl, Th) ∩ [W 1,1(Th)]3 and let Φ ∈ [C∞

0 (Ω)]3.

Ω

v · curl Φ =

T∈Th
T

v · curl Φ =

T∈Th
−
∂T

(nT × v) · Φ +

T

(curl v) · Φ

=
Ω

curlh v · Φ −

F∈Fi

h

F

[v]T · Φ (3.1) where we have used that v ∈ [W 1,1(Th)]3 so that v|T ∈ [L1(T)]3 which allows us to split the boundary integral into the integrals over the faces. Now assume that [v]T = 0 ∀F ∈ Fi

h, then

Ω

v · curl Φ =

Ω

curlh v · Φ. If we compare this to the definition of the weak derivative,

Ω

v · curl φ =

Ω

curl v · φ ∀φ ∈ [C∞

0 (Ω)]3,

we see that curl v = curlh v because Φ is arbitrary. This implies v ∈ H(curl, Ω). Conversely assume v ∈ H(curl, Ω), then relation (3.1) becomes

F∈Fi

h

F

[v]T · Φ =

Ω

curlh v · Φ −

Ω

v · curl Φ =

Ω

(curlh v − curl v) · Φ

Lemma 3.7

= Because Φ is arbitrary this implies [v]T = 0 on all F ∈ Fi

h.

Lemmas 3.5 and 3.6 are proven in a similar fashion as Lemmas 3.7 and 3.8. We refer the reader to the book of Di Pietro and Ern [6, Sect. 1.2.5] for the actual proofs.

Remark. Lemma 3.8 gives another explanation for the continuity of the tangential com-

ponent of the magnetic field H previously discussed in Sect. 2.1: In order for H to fulfil the strong form of Ampères law (2.1d) its curl must exist, hence H ∈ H(curl, Ω) (assuming that no surface currents exist). 16

SLIDE 20

3.1. Preliminaries 3.1.3. Discrete Friedrichs inequality In order to prove convergence of the DG formulation later on we need another theorem which is based on the spaces: Sh :=

ph ∈ H1

0(Ω) : ph|T ∈ Pk+1 3

(T), ∀T ∈ Th

Xh :=
vh ∈ [Pk

3(Th)]3 : (vh, grad ph) = 0, ∀ph ∈ Sh

X :=
v ∈
W 1,1(Th)

3

div v = 0 in Ω
where the space Xh contains only ”discrete divergence free” functions.

Note that if v ∈ X then (v, grad ph) = 0 for all ph ∈ Sh. Indeed, if v ∈ X then [v · nF ] = 0 over all faces F ∈ Fi

h (see [6, Lemma 1.24]) and we can write

(v, grad ph) =

T∈Th

v · grad ph =

T∈Th
F∈FT
F

(v · nT )ph =

F∈Fi

h

F

[phv] · nF = 0 where we have used that v is divergence free and that neither ph nor v jumps in normal direction. The following theorem is due to Creusé and Nicaise [5, Theorem 4.5]: Theorem 3.9. Let TH be an admissible, weakly quasi-uniform sequence of conforming meshes which consist only of tetrahedrons or hexahedrons. Then there exists a positive constant Cf > 0 independent of h such that for all h ∈ H, one has AhL2(Ω) ≤ Cf

 curlh AhL2(Ω) +

F∈Fh
F

h−1

F |[Ah]T |2

  ,

∀Ah ∈ Xh (3.2)

Remark. This is the discrete version of the first Friedrich inequality (see e.g. [13, Corol-

lary 3.51]). It is based on a discrete compactness property which has been proven for conforming/matching meshes only (see [12]). 3.1.4. Trace inequalities and polynomial approximation Lemma 3.10 (Discrete trace inequality). Let Th be an admissible mesh sequence. Then, for all h ∈ H, all vh ∈ Pk

d(Th), all T ∈ Th,

h

1/2

T vhL2(∂T) ≤ CtrN

1/2

∂

vhL2(T) (3.3) where Ctr and N∂ are both independent of h and T.

Proof. See [6, Lemma 1.46 and Eq. (1.40)]

17

SLIDE 21

3. Discontinuous Galerkin formulation

Lemma 3.11 (Polynomial approximation of curl operator). Let Th be an admissible mesh sequence and πh be the L2-orthogonal projection onto [Pk

3(Th)]3.

Then for all A ∈ [Hk+1(Th)]3, there holds curlh (A − πhA)L2(Ω) ≤ Chk|A|Hk+1(Ω) where C is independent of h.

Proof. Let u = (ux, uy, uz)T ∈ [H1(T)]3 then we have:

curl u2

L2(T) = ∂yuz − ∂zuy2 L2(T) + ∂zux − ∂xuz2 L2(T) + ∂xuy − ∂yux2 L2(T)

≤ 2|u|2

H1(T)

Applying [6, Lemma 1.58], curl (A − πhA)L2(T) ≤ C′|A − πhA|H1(T) ≤ Chk

T |A|Hk+1(T).

Summing over all T ∈ Th and using hT ≤ h for all T ∈ Th yields the assertion.

3.2. DG formulation of Eddy Current problem

In this section the symmetric-weighted-interior-penalty discretization of the eddy current IBVP (2.13) is presented and convergence to the analytical solution is proven. Many

f the lemma’s and proofs in this section are inspired by lemmas/proofs in the book of

Di Pietro and Ern [6]. In order to simplify the problem we make the additional Assumptions 3.12.

1. The computational domain Ω is a simply connected polyhedron with Lipschitz

boundary and can therefore be tiled by a mesh Th consisting of polyhedral elements T ∈ Th.

2. There exists a partition PΩ = {Ωi}1≤i≤Np of Ω such that

a) Each Ωi is a polyhedron. b) The permeability µ and the conductivity σ are constant on each Ωi and independent of time. Moreover there exist constants µmin, µmax, σmax such that 0 < µmin ≤ µ ≤ µmax and 0 < σ ≤ σmax.

3. The mesh sequence TH is compatible with the partition PΩ: For each Th ∈ TH,

each element T ∈ Th belongs to exactly one Ωi ∈ PΩ.

4. The initial condition A(t = 0) and the impressed current ji(t) are both divergence

free and the boundary traces exist for all t ∈ [0, tF ], i.e. A(t = 0), ji(t) ∈ X. 18

SLIDE 22

3.2. DG formulation of Eddy Current problem

5. For the Neumann boundary condition there holds (g × n) ∈ C1(∂Ω)2, i.e. the

tangential components of g are continuously differentiable. Assumption 2b implies that Ω2 in (2.13) is empty (only conductors are present) and therefore it suffices to solve the following IBVP: curl 1 µ curl A + σ∂A ∂t = ji

n Ω

(3.4a) A(t = 0) = A(0)

n Ω

(3.4b) 1 µ curl A × n = g × n

n ∂Ω.

(3.4c) Using the finite element space Vh := [Pk

3(Th)]3, the SWIP space semi-discretized varia-

tional form of the strong form is Find Ah ∈ C0(Vh) ∩ C1([L2(Ω)]3) such that for all t ∈ [0, tF ], all A′

h ∈ Vh,

σ∂Ah

∂t , A′

h

+ aSWIP

h

(Ah, A′

h) =

ji, A′

h

+
∂Ω

(g × n) · A′

h

(3.5a) Ah(t = 0) = πh(A(0)) (3.5b) Where g ∈ L2(∂Ω) and the SWIP bilinear form is aSWIP

h

(Ah, A′

h) =

Ω

1 µ curlh Ah · curlh A′

h −

F∈Fi

h

F

1

µ curlh Ah

ω

·

A′

h

T −
F∈Fi

h

F

1

µ curlh A′

h

ω

· [Ah]T +

F∈Fi

h

ηγµ,F hF

F

[Ah]T ·

A′

h

T

. (3.6) The last three summands of aSWIP

h

are usually referred to as consistency, symmetry and penalty terms, respectively. The penalty parameter η must be choosen big enough for convergence (see lemma 3.16). The factor γµ,F and the average weights are given by γµ,F := 2 µ1 + µ2 , ω1 := µ1 µ1 + µ2 , ω2 := µ2 µ1 + µ2 . (3.7) where µi denotes the (constant) value of µ on either side of a face F.

Remark. If we choose the conforming finite element space V ′

h such that V ′ h ⊂ H(curl, Ω),

the consistency, symmetry and penalty terms in aSWIP

h

would disappear and we would end up with the usual finite element formulation of the eddy current problem. Also notice that aSWIP

h

is a symmetric bilinear form and thus the resulting stiffness matrix is symmetric as well. 19

SLIDE 23

3. Discontinuous Galerkin formulation
Remark. The bilinear form aSWIP

h

is equivalent to the one stated by Creusé and Nicaise [5] if we set µ ≡ 1 everywhere. The main difference is, that the present form also supports discontinuous µ in a stable manner.

Remark. Because we have Neumann boundary conditions there is no need to integrate
ver the boundary faces. If we had Dirichlet boundary conditions aSWIP

h

would look the same but the face integrals would include the boundary faces. The problem (3.5) is a parabolic problem with Neumann boundary conditions. We will therefore use implicit Euler timestepping to avoid stability problems. For this we divide the time interval [0, tF ] into N equal intervals, δt = tF/N. We use the superscripts (·)(i) to denote the value of a quantity at time t(i) := i · δt so that the fully discretized version

f (3.5) reads as

Find A(i)

h ∈ Vh, i = 1, . . . , N such that for all A′ h ∈ Vh, we have

σA(i+1)

h

− A(i)

h

δt , A′

h

+ aSWIP

h

(A(i+1)

h

, A′

h) =

ji,(i+1), A′

h

+
∂Ω

(n × A′

h) · g(i+1),

(3.8a) A(0)

h

= πh(A(0)) (3.8b) 3.2.1. Convergence In this section we will prove that the approximate solution Ah obtained by Eq. (3.8) converges to the analytical solution A and give an upper bound for the error. Of course A must exist and we additionally require Assumption 3.13. We assume that A ∈ C0(V∗)∩C2(L2(Ω)3) where V∗ := H(curl, Ω)∩ H2(PΩ) and that A is a solution of the strong form (3.4). The first assumption is needed to assure that the traces of A on the element faces exist which in turn allows us to extend the meaning of aSWIP

h

to the space V ∗: aSWIP

h

: V ∗

h × Vh → R where V ∗ h := Vh ∪ V ∗. (The usual trace theorem for H(curl, T) gives a

trace in H−1/2

t

(∂T) which doesn’t suffice: We need the trace to be in at least L2(∂T) so that the surface integral can be split into face integrals, cf. [13, Remark 3.32] and [6, Remark 1.25]) Together with the second assumption this allows us to prove consistency: Lemma 3.14 (Consistency). The bilinearform aSWIP

h

is consistent: If A ∈ C0(V ∗

h ) ∩

C1(L2(Ω)) is the strong solution of the IBVP (3.4) then aSWIP

h

(A(t), A′

h) =

Ω

curl 1 µ curl A · A′

h +

∂Ω

(g × n) · A′

h.

(3.9) 20

SLIDE 24

3.2. DG formulation of Eddy Current problem Moreover,

σ∂A(t)

∂t , A′

h

+ aSWIP

h

(A(t), A′

h) =

ji, A′

h

+
∂Ω

(g × n) · A′

h

(3.10) for all A′

h ∈ Vh, all t ∈ [0, tF ], all h ∈ H.

Proof. Because A(t) ∈ H(curl, Ω), all tangential jumps in the definition of aSWIP

h

are zero and therefore the symmetry and the penalty term drop out immediately and we are left with: aSWIP

h

(A, A′

h) =

Ω

1 µ curl A · curlh A′

h −

F∈Fi

h

F

1

µ curl A

ω

·

A′

h

T

(3.11) Next we need the identity [a × b] · n = (a1 × b1 − a2 × b2) · n, = ((ω1a1 + ω2a2) × (b1 − b2) + (a1 − a2) × (ω1b1 + ω2b2)) · n, = − [b]T · {a}ω + [a]T · {b}ω , (3.12) where {b}ω := (ω2b1 + ω1b2). Note that this inequality is valid for any weights ω1, ω2.

Eq. (3.11) can be written as:

aSWIP

h

(A, A′

h) =

Ω

1 µ curl A · curlh A′

h +

F∈Fi

h

F

1

µ curl A × A′

h

· nF −
F∈Fi

h

F

1

µ curl A

T

·

A′

h

ω

The last term in this equation is zero because A is a strong solution of IBVP (3.4) and therefore 1

µ curl A ∈ H(curl, Ω). The integral over all faces can be split into an integral

ver all elements,

aSWIP

h

(A, A′

h) =

T∈Th
T

1 µ curl A · curl A′

h +

T∈Th
∂T

1

µ curl A × A′

h

· nT −
∂Ω

1

µ curl A × A′

h

· n

=

T∈Th
T

curl 1 µ curl A · A′

h −

∂Ω
n × 1

µ curl A

· A′

h.

Where we have used an integration by parts formula for H(curl, Ω) functions in the last step (see e.g. [13, Theorem 3.29]). From this Eq. (3.9) follows immediately. Eq. (3.10) is obtained by multiplying Eq. (3.4a) by A′

h, integrating over Ω and using Eq. (3.9).

21

SLIDE 25

3. Discontinuous Galerkin formulation

In the following we will also prove boundedness and coercivity for the discrete bilinear form aSWIP

h

. With the help of the Lax-Milgram lemma this would allow us to prove the well posedness of the time independent variational problem arising from (3.5) by setting σ ≡ 0. However we want to prove a convergence result for the time dependent problem so further steps are needed. But let us first define the (semi-)norms with respect to which we will prove coercivity and boundedness: |A|SWIP :=

µ−1/2 curlh A
2

L2(Ω) + |A|2 J,µ

1/2

, |A|J,µ : =

 

F∈Fh

γµ,F hF [A]T 2

L2(F)

 

1/2

, |A|SWIP,* :=

 |A|2

SWIP +

T∈Th

hT

µ−1/2 curlh A
T
2

L2(∂T)

 

1/2

.

Remark. The semi-norm |·|SWIP and the jump semi-norm |·|J,µ are well defined for all

functions ∈ H(curl, Ω). This is not true for the last semi norm, |·|SWIP,*, because it requires us to give a meaning to trace of curlh A on ∂T. We avoid this problem by assuming A ∈ V ∗ (Assumption 3.13) which allows for a suitable definition of the trace. Note that the defined semi-norms are not equivalent to L2(Ω) norms if A ∈ V ∗. However because the evolution equation (3.8a) only acts on the discrete divergence free components of Ah we will be able to prove convergence in the L2(Ω) norm. (A mathematically more rigorous statement follows later.) Before we prove coercivity and boundedness we state a useful lemma to bound the consistency and symmetry terms later on: Lemma 3.15. For all (Ah, A′

h) ∈ V∗ h × V∗ h there holds:

F∈Fi

h

F
µ−1 curlh Ah
ω ·

A′

h

T
≤

 

T∈Th

hT 2

µ−1/2 curlh Ah
T
2

L2(∂T)

 

1/2

A′

h

J,µ
Proof. By the Cauchy Schwarz inequality,
F
µ−1 curlh Ah
ω ·

A′

h

T
≤
F
ω1

µ1 curlh Ah,1 + ω2 µ2 curlh Ah,2

21/2
A′

h

T
2

L2(F)

1/2

22

SLIDE 26

3.2. DG formulation of Eddy Current problem

   

⇒ Owing to the Cauchy Schwarz inequality,

aibi
2 ≤
a2

i

b2

i

ω1

µ1 curlh Ah,1 + ω2 µ2 curlh Ah,2

2

≤

ω2

1

µ1 + ω2

2

µ2

1

µ1 (curlh Ah,1)2 + 1 µ2 (curlh Ah,2)2



  

≤

µ−1/2

1

curlh Ah,1

2

L2(F) +

µ−1/2

2

curlh Ah,2

2

L2(F)

1/2

ω2

1

µ1 + ω2

2

µ2

A′

h

T
2

L2(F)

1/2     

⇒ Now use that, ω2

1

µ1 + ω2

2

µ2 = µ1 (µ1 + µ2)2 + µ2 (µ1 + µ2)2 = 1 µ1 + µ2 = 1 2γµ,F , and sum over all edges

    

F∈Fi

h

F
µ−1 curlh Ah
ω ·

A′

h

T
≤

  

F∈Fi

h

γµ,F hF

A′

h

T
2

L2(F)

  

1/2

  

F∈Fi

h

hF 2

µ−1/2

1

curlh Ah,1

2

L2(F) + hF

2

µ−1/2

2

curlh Ah,2

2

L2(F)

  

1/2

≤

 

T∈Th

F∈FT

hF 2

µ−1/2 curlh Ah
T
2

L2(F)

 

1/2

A′

h

J,µ .
bserving that hF ≤ hT yields the assertion.

Lemma 3.16 (Discrete Coercivity). The bilinear form aSWIP

h

is coercive w.r.t. to the |·|SWIP semi-norm: For all η > C2

trN∂/

√ 2, all Ah ∈ Vh there exists a constant Cstab > 0

such that for all h ∈ H aSWIP

h

(Ah, Ah) ≥ Cstab |Ah|2

SWIP

with Cstab =

η−C2

trN∂/ √ 2

1−η

. The constant Ctr and N∂ stem from lemma 3.10 and are independent of h.

Proof. Lemma 3.15 gives us an estimate on the consistency and symmetry term of aSWIP

h

. The factor |Ah|J,µ appears in the definition of the |·|SWIP norm hence only an estimate 23

SLIDE 27

3. Discontinuous Galerkin formulation

for the other factor on the right-hand side of lemma 3.15 is needed. Owing to the trace inequality 3.10 (Ah ∈ Vh),

T∈Th

hT 2

µ−1/2 curlh Ah
T
2

L2(∂T) ≤

T∈Th

C2

trN∂

2

µ−1/2 curlh Ah
2

L2(T)

= C2

trN∂

2

µ−1/2 curlh Ah
2

L2(Ω) .

(3.13) Hence, aSWIP

h

(Ah, Ah) ≥

µ−1/2 curlh Ah
2

L2(Ω) − Ctr

2N∂
µ−1/2 curlh Ah
L2(Ω) |Ah|Jµ

+

F∈T i

h

ηγµ,F hF [Ah]T 2

L2(F) .

Next we consider the inequality x2−2βxy+ηy2 ≥ η−β2

1−η (x2+y2) which holds for arbitrary

β, η, x, y,. (It can be derived from from (βx − ηy)2 + (x − βy)2 ≥ 0.) Applying this to the above estimate, aSWIP

h

(Ah, Ah) ≥ η − C2

trN∂/

√ 2

1 − η

µ−1/2 curlh Ah
2

L2(Ω) + |Ah|2 J,µ

= Cstab |Ah|SWIP .

Note that Cstab > 0 only if the penalty parameter η > C2

trN∂/

√ 2 which completes the

proof. Lemma 3.17 (Boundedness). There exists a constant Cbnd independent of h such that for all (Ah, A′

h) ∈ V ∗ h × Vh, all h ∈ H

aSWIP

h

(Ah, A′

h) ≤ Cbnd |Ah|SWIP,*

A′

h

SWIP
Proof. Recall the definition of the bilinear form,

aSWIP

h

(Ah, A′

h) =

Ω

1 µ curlh Ah · curlh A′

h −

F∈Fi

h

F

1

µ curlh Ah

ω

·

A′

h

T −
F∈Fi

h

F

1

µ curlh A′

h

ω

· [Ah]T +

F∈Fi

h

ηγµ,F hF

F

[Ah]T ·

A′

h

T

= T1 + T2 + T3 + T4 24

SLIDE 28

3.2. DG formulation of Eddy Current problem Next we need to bound the individual terms by the respective norms: |T1| ≤

µ−1/2 curlh Ah
L2(Ω)
µ−1/2 curlh A′

h

L2(Ω) ≤ |Ah|SWIP
A′

h

SWIP ,

|T2|

Lemma 3.15

≤ |Ah|SWIP,*

A′

h

SWIP ,

|T3|

Lemma 3.15

≤

 

T∈Th

hT 2

µ−1/2 curlh A′

h

T
2

L2(∂T)

 

1/2

|Ah|J,µ ,

Eq. (3.13)

≤

C2

trN∂

2

µ−1/2 curlh A′

h

L2(Ω) |Ah|J,µ ≤ C
A′

h

SWIP |Ah|SWIP

|T4| ≤

F∈Fi

h

ηγµ,F hF [Ah]T L2(F)

A′

h

T
L2(F) ≤ η |Ah|J,µ
A′

h

J,µ

≤ η |Ah|SWIP

A′

h

SWIP ,

where we have used the Cauchy-Schwarz inequality multiple times to bound T1 and T4. Thus we have proved the coercivity and boundedness of aSWIP

h

. Consider for a moment the discrete problem which is obtained from Eq. (3.5a) by setting σ ≡ 0 in Ω. We could now use the Lax-Milgram Lemma to prove the well posedness of this problem if we add suitable boundary conditions and require Ah to be discrete divergence free. The last condition is needed in order to render Ah unique and implies that the |·|SWIP and |·|SWIP,* semi-norms are actually full norms on the space Vh. This is a more rigorous explanation for the introduction of the gauge condition which was already discussed in

Chap. 2.

However, our goal is to proof convergence of the discrete evolution equation (3.8a) as h → 0. Therefore an expression for the evolution of the approximation error A(i) − A(i)

h

is needed. For this we split the total error, A(i)

h − A(i) = A(i) h − πhA(i)

ξ(i)

h

− (A(i) − πhA(i))

ξ(i)

π

. The projection error ξ(i)

π

can always be estimated using polynomial approximation re-

sults. Hence only the development of the term ξ(i)

h

is of importance for the moment. Lemma 3.18 (Evolution of ξ(i)

h ). Assume the exact solution A ∈ C0(V ∗ h ) ∩ C2(L2(Ω)),

then

σξ(i+1)

h

− ξ(i)

h

δt , A′

h

= −aSWIP

h

(ξ(i+1)

h

, A′

h) + α(i+1)

∀A′

h ∈ Vh,

25

SLIDE 29

3. Discontinuous Galerkin formulation

where α(i+1) := aSWIP

h

(ξ(i+1)

π

, A′

h) −

σπhθ(i+1), A′

h

and θ(i+1) := δt−1 t(i+1)

t(i)

(t(i) − t)∂2A(t)

∂t2

dt

Proof. Owing to Taylor’s theorem,

A(i) = A(i+1) − δt∂A(i+1) ∂t −

t(i+1)

t(i)

(t(i) − t)∂2A(t) ∂t2 dt. Projecting this equation into Vh we get πhA(i) = πhA(i+1) − δtπh ∂A(i+1) ∂t − δtπhθ(i+1), Now multiply on both sides with σ and test with A′

h ∈ Vh,

σπhA(i), A′

h

=
σπhA(i+1), A′

h

− δt
σπh

∂A(i+1) ∂t , A′

h

− δt
σπhθ(i+1), A′

h

,

and owing to the consistency of aSWIP

h

,

σπh

A(i+1) − A(i) δt , A′

h

= −aSWIP

h

(A(i+1), A′

h) +

ji,(i+1), A′

h

+
∂Ω

(g(i+1) × n) · A′

h +

σπhθ(i+1), A′

h

Finally we subtract the above equation from the discrete evolution equation (3.8a) to

get

σξ(i+1)

h

− ξ(i)

h

δt , A′

h

= −aSWIP

h

A(i+1)

h

− A(i+1), A′

h

−
σπhθ(i+1), A′

h

= −aSWIP

h

ξ(i+1)

h

, A′

h

+ α(i+1)

Next we are going to prove that the solution A(i)

h

stays discrete divergence free if the initial value A(0) is divergence free(cf. Sec. 3.1.3), Lemma 3.19. All solutions A(i+1)

h

f the evolution equation (3.8) are discrete divergence

free, i.e. they belong to the space Xh. Moreover we have for the exact solution A(t) of the IBVP (3.4) (A(t), grad p) = 0 for all p ∈ Sh, t ∈ [0, tF ]. 26

SLIDE 30

3.2. DG formulation of Eddy Current problem

Proof. Observe that if p ∈ Sh, then grad p ∈ H(curl, Ω). Indeed, p can be viewed

as a discrete 0-form and grad p is therefore a discrete one form and thus tangentially

continuous. See [3] for a detailed explanation. Hence we have [ grad p]T = 0 over all

Faces F ∈ Fi

h and grad p = grad hp. Therefore the following terms cancel out:

aSWIP

h

(A(i+1)

h

, grad p) = 0 by definition of aSWIP

h

ji,(i+1), grad p
= 0

by assumption 3.12

∂Ω

(g × n) · grad p =

∂Ω

(n × ( grad p − grad p · n)) · g = −

∂Ω

g · curlΓp = −

∂Ω

p curlΓ g = 0 where curlΓ denotes the vectorial surface curl and curlΓ the scalar surface curl (cf. [13,

Sect. 3.4]). In the last step we applied a variant of Stokes theorem (see e.g. [13, Corollary

3.21]) and used the fact that p = 0 on ∂Ω. Let us now substitute A′

h = grad p into Eq. (3.8a),

σA(i+1)

h

, grad p

=
σA(i)

h , grad p

.

Summing the above equation for i = 0, . . . N−1 we see that

A(N)

h

, grad p

=
A(0)

h , grad p

=

(πhA(t = 0), grad p) = (A(t = 0), grad p) = 0 which proves the first assertion. For the second claim we multiply Eq. (3.4a) by grad p (p ∈ Sh) and split the integral over all elements so we can use the partial integration formula (see e.g. [13, Theorem 3.29]),

curl 1

µ curl A(t), grad p

+
σ∂A(t)

∂t , grad p

=
ji, grad p
T∈Th
T

curl 1 µ curl A(t) · grad p + ∂(σA(t), grad p) ∂t = 0 −

T∈Th
F∈FT
F

1

µ curl A(t) × nT

· grad p
T

+ ∂(σA(t), grad p) ∂t = 0. Now applying the relation (3.12) and noting that

µ−1 curl A(t)

T = 0 and [ grad p]T =

∂(σA(t), grad p) ∂t = 0. Observe that A(0) ∈ X and therefore (A(t), σ grad p) = 0 for all p ∈ Sh which yields the assertion. 27

SLIDE 31

3. Discontinuous Galerkin formulation

Lemma 3.19 is important because it implies that ξ(i+1)

h

∈ Xh as well. Indeed we have

ξ(i+1)

h

, grad p

=
A(i+1)

h

− πhA(i+1), grad p

= 0 +
A(i+1), grad p
= 0 for all

p ∈ Sh. Lemma 3.20 (Energy estimate). Under the assumptions of lemma 3.18 assume in addition that the mesh sequence TH is admissible, weakly quasi-uniform and matching, then there holds for i = 0, . . . , N − 1,

√σξ(i+1)

h

2

L2(Ω) +

√σ(ξ(i+1)

h

− ξ(i)

h )

2

L2(Ω) + Cstabδt

ξ(i+1)

h

2

SWIP ≤

√σξ(i)

h

2

L2(Ω)

+ Cδt

ξ(i+1)

π

2

SWIP,* + C2 Aδt2

,

where CA = maxt∈[t(i),t(i+1)]

∂2A(t)

∂t2

L2(Ω).
Proof. Choose A′

h = ξ(i+1) h

in lemma 3.18 and multiply with δt,

√σξ(i+1)

h

2

L2(Ω) −

σξ(i)

h , ξ(i+1) h

= −δtaSWIP

h

(ξ(i+1)

h

, ξ(i+1)

h

) + δtaSWIP

h

(ξ(i+1)

π

, ξ(i+1)

h

− δt

σπhθ(i+1), ξ(i+1)

h

,

Now apply the relation −2ab = (a − b)2 − a2 − b2 to the term

σξ(i)

h , ξ(i+1) h

,
√σξ(i+1)

h

2

L2(Ω) −

√σξ(i)

h

2

L2(Ω) +

√σ(ξ(i+1)

h

− ξ(i)

h )

L2(Ω) = −2δtaSWIP

h

(ξ(i+1)

h

, ξ(i+1)

h

)+ 2δtaSWIP

h

(ξ(i+1)

π

, ξ(i+1)

h

) − 2δt

σθ(i+1), ξ(i+1)

h

,

use the coercivity and boundedness of aSWIP

h

as well as the Cauchy Schwarz inequality,

√σξ(i+1)

h

2

L2(Ω) +

√σ(ξ(i+1)

h

− ξ(i)

h )

L2(Ω) + 2δtCstab
ξ(i+1)

h

2

SWIP ≤

√σξ(i)

h

2

L2(Ω) +

2Cbndδt

ξ(i+1)

π

SWIP,*
ξ(i+1)

h

SWIP + 2δt
σθ(i+1)
L2(Ω)
ξ(i+1)

h

L2(Ω) .

Because ξ(i+1)

h

∈ Xh we can use the discrete Friedrich’s inequality 3.2 for the last term,

√σξ(i+1)

h

2

L2(Ω) +

√σ(ξ(i+1)

h

− ξ(i)

h )

L2(Ω) + 2δtCstab
ξ(i+1)

h

2

SWIP ≤

√σξ(i)

h

2

L2(Ω) +

2δt

ξ(i+1)

h

SWIP ·
Cbnd
ξ(i+1)

π

SWIP,* + C′

F

σθ(i+1)
L2(Ω)
28

SLIDE 32

3.2. DG formulation of Eddy Current problem using Youngs inequality, 2ab ≤ a2/ε + εb2, with ε = C−1

stab we finally get

√σξ(i+1)

h

2

L2(Ω) +

√σ(ξ(i+1)

h

− ξ(i)

h )

L2(Ω) + δtCstab
ξ(i+1)

h

2

SWIP ≤

√σξ(i)

h

2

L2(Ω) +

δtCstab

Cbnd
ξ(i+1)

π

SWIP,* + C′

F

σθ(i+1)
L2(Ω)

2

. The only thing left to complete the proof is to give an upper bound on

σθ(i+1)
L2(Ω),
σθ(i+1)
2

L2(Ω) ≤ σ2 maxδt−2

Ω

t(i+1)

t(i)

(t(i) − t)∂2A(t) ∂t2 dt

2

dx, ≤ σ2

maxδt−2

Ω

t(i+1)

t(i)

(t(i+1) − t)2dt

  t(i+1)

t(i)

∂2A(t)

∂t2

2

dt

  dx,

≤ 1 3σ2

maxδt

 δt

max

t∈[t(i),t(i+1)]

∂2A(t)

∂t2

2

L2(Ω)

  ≤ C′

Aδt2,

where C′

A depends on the analytical solution.

Lemma 3.20 is essential in proving the convergence of the numerical solution Ah to the exact solution as h → 0. We are now finally in the position to state the main result of this chapter. Theorem 3.21 (Convergence). Let TH be an admissible, weakly quasi-uniform, matching mesh sequence. Assume assumptions 3.12 and 3.13 hold, A ∈ C0(Hk+1(PΩ)), Ah ∈ [Pk

3(TH)]3 then for all h ∈ H , all integer k > 0,

√σ(A(N) − A(N)

h

)

L2(Ω) +
Cstabδt

N

i=1
A(i) − A(i)

h

2

SWIP

1/2

≤ Ct

1/2

F

C1hk + C2δt
where C1 = maxt∈[0,tF ] |A(t)|Hk+1(Ω) and C2 = maxt∈[0,tF ]
∂2A(t)

∂t2

L2(Ω). The constants

C1,C2 and C are independent of h and δt.

Proof. In the following C will denote an arbitrary, positive constant which is not always

the same. Let us first proof a polynomial approximation result to bound the term

ξ(i)

π

SWIP,* in terms of h,
ξ(i)

π

2

SWIP,* =

µ− 1

2 curlh ξ(i)

π

2

L2(Ω) +

F∈Fh

γµ,F hF

ξ(i)

π

T
2

L2(F) +

T∈Th

hT

µ− 1

2 curlh ξ(i)

π

T
2

L2(∂T)

= T1 + T2 + T3 29

SLIDE 33

3. Discontinuous Galerkin formulation

Owing to lemma 3.11 the first term can easily be bounded by T1 ≤ Ch2k

A(i)
2

Hk+1(Ω).

For the second term we have T2 ≤

T∈Th
F∈FT

2γµ,F hF

ξ(i+1)

π

T
2

L2(F) ≤ Ch2k T

A(i)
2

Hk+1(Ω)

where we have used [6, Lemma 1.59] in the last step. Applying the same lemma to the third term T3, T3 ≤ C

T∈Th
F∈FT

hT

3

j=1
grad (ξ(i)

π )j

T
2

L2(Ω) ≤ Ch2k

A(i)
2

Hk+1(Ω)

Hence

ξ(i)

π

SWIP,* ≤ Capphk
A(i)
Hk+1(Ω). Summing up the energy estimate from lemma

3.20 over i = 0, . . . , N − 1 and dropping the non-negative term

√σ(ξ(i+1)

h

− ξ(i)

h )

2

L2(Ω),

√σξ(N)

h

2

L2(Ω) + Cstabδt N

i=1
ξ(i)

h

2

SWIP ≤ C

C2

1h2k + C2 2δt2

.

      

⇒ Owing to the triangle inequality,

√σ(A(N) − A(N)

h

)

2

L2(Ω) ≤ 2

√σξ(N)

h

2

L2(Ω) + 2

√σξ(N)

π

2

L2(Ω) N

i=1
A(i) − A(i)

h

2

SWIP ≤ 2 N

i=1
ξ(i)

h

2

SWIP + 2 N

i=1
ξ(i)

π

2

SWIP

      

1 2

√σ(A(N) − A(N)

h

)

2

L2(Ω) + 1

2Cstabδt

N

i=1
A(i) − A(i)

h

2

SWIP ≤ C

C2

1h2k + C2 2δt2

+

√σξ(N)

π

2

L2(Ω) + δt N

i=1
ξ(i)

π

2

SWIP .

The two last terms are also bounded in terms of h2k, δt

N

i=1
ξ(i)

π

2

SWIP ≤ N

i=1
ξ(i)

π

2

SWIP,* ≤ C · C2 1h2k,

√σξ(N)

π

2

L2(Ω) ≤ √σC′ apph2k+2

A(N)
Hk+1(Ω),

where the last inequality is a consequence of [6, Lemma 1.58]. This completes the proof. 30

SLIDE 34

3.2. DG formulation of Eddy Current problem The convergence theorem 3.21 is optimal w.r.t. the |·|SWIP norm, i.e. the optimal convergence rate k is achieved. The only flaw is that for k = 1 the H2 semi-norm appears on the right-hand side instead of a term of the form

1

µ curlh A

H1(Ω). However

this is not limitation in light of assumption 3.13. In contrast to the SWIP-norm the bound on the L2 norm is suboptimal by one power in h. However it should be possible to get a sharper bound by following the idea of Wheeler (cf. [6, Sect. 4.7.5]). It is possible to derive the same error estimate if A is approximated with edge functions and not general polynomials. In this case one can use Theorem 5.41 from [13] in the proof of Thm. 3.21. The biggest flaw, with respect to this thesis, is that Thm. 3.21 is only valid for matching meshes, i.e. hanging nodes are not allowed. This stems from the discrete compactness property which is only valid for matching meshes (see [12]) and on which the discrete Friedrich inequality (Thm. 3.9) depends. A convergence theorem similar to Thm. 3.21 can be derived for dirichlet boundary conditions by following a similar approach. The only difference in the aSWIP

h

bilinear form is that the sums over all faces also include boundary faces. Moreover it should be possible to relax the condition on the Neumann boundary condition from g×n ∈ C1(∂Ω) to g × n ∈ H1(∂Ω)2 by splitting the surface integral in the proof of lemma 3.15 into integrals over the faces of each element. Finally we want to stress that this convergence result is only valid for h-convergence. For p refinement further investigations are needed. 31

SLIDE 35

4. Aspects of the Implementation

So far we have always looked at the 3D eddy current problem. In chapter 2 we have derived the model from a physical point of view and in chapter 3 we have given a proof

f convergence under the assumptions 3.12 and 3.13. We are now going to simplify the

problem from chapter 3 by

1. reducing it to two dimensions (Sec. 4.1), and by
2. assuming σ > 0 everywhere i.e. the whole computational domain is filled by a

conductor. The second assumption is not a real restriction: The electromagnetic fields in free space can be calculated after the fields inside the conductors have been calculated (cf. [2]). On the other hand we will extend the applicability of the model by dropping the assumption of Ω being a simply connected polyhedron. When we discuss the numerical results in chapter 5 we will consider settings where Ω or a subdomain Ωi ∈ PΩ has the shape of a circle and therefore no polyhedral mesh can exist. To handle such geometries in a consistent way we introduce special elements with circle shaped boundaries (Sec. 4.2). This chapter starts off by presenting four different ways to discretize the 3D Eddy current problem from the previous chapter under the assumption that the whole system is at

rest. Following this we present some implementation details. Finally we merge the 2D

variational formulation obtained from Sec. 4.1 with the transformation laws of Sec. 2.4 in Sec. 4.3.

4.1. Reduction to 2D

Our goal is to reduce the temporal gauged eddy current formulation 2.13 and the H- formulation into a meaningful two-dimensional representation. There are principally two different ways of doing this: The transverse electric (TE) and the transverse magnetic (TM) formulation (cf. [9, Sec. 8.2], [4]). In both cases one assumes that the 2D domain Ω lies in the x − y plane and that it is the cross-section of a 3D unbounded cylinder which is aligned with the z-axis. Further the material coefficients µ and σ as well as the

SLIDE 36

4.1. Reduction to 2D source current ji are independent of the z-component. In the TE model it is assumed that the impressed current ji lies completely in the x − y plane, ji = (ji

x, ji y, 0). In contrast the TM model assumes that ji is aligned with the

z-axis, i.e. ji = (0, 0, ji

z). In both cases the form of the remaining electromagnetic fields

is easily derived. We have summarized them in table 4.1. Field quantity TM-formulation TE-formulation ji resp. jf (0, 0, ji

z), (0, 0, jf z )

(ji

x, ji y, 0), (jf x, jf y , 0)

E (0, 0, Ez) (Ex, Ey, 0) B resp. H (Bx, By, 0), (Hx, Hy, 0) (0, 0, Bz), (0, 0, Hz) A resp. A∗ (0, 0, Az), (0, 0, A∗

z)

(Ax, Ay, 0), (A∗

x, A∗ y, 0)

Table 4.1.: The 2D equivalents of most electromagnetic quantities for the TE and TM formulation. In order to state the 2D variational form of the eddy current models it is useful to introduce the following 2D scalar and vector curl operators: curl2D A := ∂Ay ∂x − ∂Ax ∂y , (2D Scalar curl) curl2D Az :=

∂Az

∂y

−∂Az

∂x

.

(2D Vector curl) Let us further introduce the 2D transverse respectively longitudinal version of the symmtric- weighted-interior-penalty bilinearform (3.6):

µaSWIP h,T

(Ah, A′

h) :=

Ω

1 µ curl2D,h Ah · curl2D,h A′

h −

F∈Fi

h

F

1

µ curl2D,h Ah

ω

·

A′

h

T

−

F∈Fi

h

F

1

µ curl2D,h A′

h

ω

· [Ah]T +

F∈Fi

h

ηγµ,F hF

F

[Ah]T ·

A′

h

T

µaSWIP h,L

(Az, A′

z) :=

Ω

1 µ grad hAh,z · grad hA′

h,z −

F∈Fi

h

F

1

µ grad hAh,z

ω

· nF

A′

h,z

−
F∈Fi

h

F

1

µ grad hA′

h,z

ω

· nF [Ah,z] +

F∈Fi

h

ηγµ,F hF

F

[Ah,z]

A′

h,z

where the tangential jump for the µaSWIP

h,T

bilinearform is defined by [Ah]T = nF,xAy − nF,yAx. The broken, scalar curl operator curl2D,h is defined analogous to the broken gradient in Def. 3.3. Let us further define the bilinearforms µaSWIP

h,T,0 respectively µaSWIP h,L,0

33

SLIDE 37

4. Aspects of the Implementation

in the same way but such that the sums of the consistency, symmetry and penalty terms include the boundary faces of the mesh Th. We need those extended forms in the next section for solving IBVP with homogeneous dirichlet boundary conditions. Note that the longitudinal bilinearform µaSWIP

h,L,0 is exactly the same as the one given by

Di Pietro and Ern [6, Eq. (4.64)] for the Poisson problem with heterogeneous diffusion. This agrees with the fact that the H-model in the TE formulation is equivalent to the 2D heat equation.

Tranverse Electric formulation (TE)

By having a peek at the variational formulation 3.8 it is easy to state the 2D TE version for the temporal gauged potential formulation (2.13): Find A(j)

h

∈ Rk(Th), j = 1, . . . , N, such that for all A′

h ∈ Rk, all i = 0, . . . , N − 1,

σA(i+1)

h

− A(i)

h

δt , A′

h

+µ aSWIP

h,T

(A(i+1)

h

, A′

h) =

ji,(i+1), A′

h

+
∂Ω

(n × A′

h) · g(i+1)

(4.1a) A(0)

h

= πh(A(t = 0)) (4.1b) where Rk(Th) is the space of Nedelec elemens of the first kind of order k over the mesh Th (see [13, Sec. 5.5]).

Remark. Instead of the space Rk we could have chosen [Pk

d]2. However in this thesis

we want to use the DG face integrals only along the interface between two bodies. Everywhere else we use H(curl, Ω) conforming basis functions and thus edge elements are the natural choice.

Remark. The variational formulation 4.1 is only correct for systems at rest. We will

state the correct variational formulation for moving systems in Sec. 4.3. In a similar way we can readily state the variational formulation for the H-model (2.16) with homogeneous dirichlet boundary conditions (g = 0) in TE-form: Find H(j)

h,z ∈ Pk 1(Th), j = 1, . . . , N, such that for all H′ h,z ∈ Pk 1(Th), all i = 0, . . . , N − 1,

 µ

H(i+1)

h,z

− H(i)

h,z

δt , H′

h,z

  +σ aSWIP

h,L,0 (H(i+1) h,z

, H′

h,z) =

curl2D,h

1 σji,(i+1), H′

h,z

,

(4.2a) H(0)

h,z = πh(Hz(t = 0)).

(4.2b) Note how the Neumann source term on the right-hand side is missing. Inhomogeneous dirichlet boundary conditions could easily be implemented by adding another term on the right-hand side. The idea is that the jumps and averages on the boundary faces are redefined using the dirichlet data (use symmetric averages with ω1 = ω2 = 1/2). 34

SLIDE 38

4.2. Details of the imlementation.

Transverse Magnetic formulation (TM)

The variational formulation of the temporal gauged potential model as well as the H model are easily derived in the same manner as for the TM-formulation. For the sake of completeness we will state them here. For the temporal gauged potential model (2.13) we have: Find A(j)

h,z ∈ Pk 1(Th), j = 1, . . . , N such that for all A′ h,z ∈ Pk 1(Th), all i = 0, . . . , N − 1,

 σ

A(i+1)

h,z

− A(i)

h,z

δt , A′

h

  +µ aSWIP

h,L

(A(i+1)

h,z

, A′

h,z) =

ji,(i+1)

z

, A′

h,z

+
∂Ω

g(i+1)A′

h,z,

(4.3a) A(0)

h,z = πh(Az(t = 0)).

(4.3b) Similarly the TM-version of the H-formulation (2.16) is, Find H(j)

h

∈ Rk(Th), j = 1, . . . , N, such that for all H′

h ∈ Rk(Th), all i = 0, . . . , N − 1,

µH(i+1)

h

− H(i)

h

δt , H′

h

+σ aSWIP

h,T,0 (H(i+1) h

, H′

h) =

curl2D,h

1 σji,(i+1)

z

, H′

h

,

(4.4a) H(0)

h

= πh(H(t = 0)). (4.4b)

4.2. Details of the imlementation.

In the previous section we have stated four 2D variational formulations. In the next sections we are discussing some details of their implementation on the Computer. In particular we will focus on the non-conforming mesh data structure. For a (almost complete) documentation of the code we refer the reader to appendix B and the Doxygen documentation.

General remarks

All the code has been written in C++11 and is based on the Fi- nite Element library NGSolve (http://sourceforge.net/projects/ngsolve/). The latter builds in turn on the mesh generator Netgen which also provides the graphical interface for meshing and displaying the solution. Most of the code which is contained in the NonConforming-library has been tested for correctness using the unit testing framework CUTE developed by Sommerlad and Graf [17]. 35

SLIDE 39

4. Aspects of the Implementation

One of the big problems which we faced during the development of the code was the lack

f support for nonconforming meshes, both by NGSolve as well as Netgen. Although

it’s possible to store nonconforming meshes in the data structures of Netgen, NGSolve doesn’t support this in any way (in particular the Discontinuous Galerkin components in NGSolve). For this reason and some usability reasons the Author decided to come up with a complete remake of the mesh data structure. This data structure is discussed in detail in the next subsection.

Remark. Because the code builds on NGSolve a lot of it’s functionality was reused. In

particular all the shape functions and most of the element mappings stem from NGSolve. The matrix/vector assembly had to be reimplemented in order to support the nonconforming mesh data structure (this was done in the classes NC_Bilinearform resp. NCLinearForm). 4.2.1. The mesh data structure The newly introduced NonConformingMesh class is in a way a different representation of the data which is available through the MeshAccess interface of NGSolve. The idea is that both data structures coexist at the same time and that they are syn-

chronized. The only difference between the two is the representation of the data: The

NonConformingMesh class provides the user with additional information about hanging nodes, elements which are not matching etc. This information could be computed directly from the MeshAccess class but in a rather complicated way. In other words the NonConformingMesh class can be seen as a buffer for storing additional information about the nonconforming nature of the mesh. Internally the NonConformingMesh stores references to a number of other classes which represent the actual geometric objects. We will refer to the entirety of all those classes as the mesh data structure from now on. Figure 4.1 shows all the involved classes. The mesh data structure has been constructed such that it can be extended to three dimensions easily. However the price for this is that the model is more complicated than what is needed for a purely 2D setting. Let us start with the geometric concepts of surfaces, curves and points (we will provide a concrete example a bit later. For now it is a good idea to have Fig. 4.1 at hand while following the description). These are all implemented in the classes Surface<TRIG>, Curve<LINE>, Curve<CIRC_SEGM> and Point. E.g. the class Surface<TRIG> represents a triangle in 3D/2D space. Notice that these geometric classes derive from the purely abstract base classes (interfaces) ISurface, ICurve. The interfaces define a common set of methods related to the geometric nature of surfaces or curves. This is essential because the matrix assembling routines presented later make no difference between circle segments or lines, they treat them in exactly the same way (loose coupling). 36

SLIDE 40

4.2. Details of the imlementation. The Point class doesn’t have such an interface because a point in space has always the same geometric properties. A second layer of mesh elements has been built on top of these purely geometric classes. At the top level a 2D NonConformingMesh consists of an arbitrary number of IFace

bjects. The present implementation supports triangles through the class Face<TRIG>

which implements IFace. Note that Face<TRIG> is essentially a triangle (implementation is inherited from Surface<TRIG>!) with additional properties. Most notably an IFace object contains references to IFaceSegment objects. In 2D each IFace object ”contains” exactly one IFaceSegment object and a one-to-one relationship could be established between them. In 3D this is not the case as should become clearer later on. In a similar fashion each ISurface object contains references to IEdge objects. The union of all those IEdge objects represents the geometric boundary of the ISurface

bject. E.g. each triangle references exactly three Edge<LINE> objects. However each

IEdge object belongs to exactly one IFace and/or IFaceSegment object. In fact in 2D every IEdge object corresponds to exactly one face and its facesegment. An edge cannot be associated with two different triangles! In order to connect the different faces of the mesh with one another IEdgeSegments have been introduced: In 2D each IEdgeSegment belongs to exactly one IEdge object if it lies on ∂Ω. Otherwise it belongs to exactly two IEdge objects. Vice versa each IEdge object consists of an arbitrary number of IEdgeSegment objects. This is required for non-conforming meshes if, for instance, the two sides of a triangle do not

match. If the mesh is matching, each IEdge object contains exactly one IEdgeSegment
bject but each IEdgeSegment can correspond to one (boundary edge) or two (inner

edge) IEdge objects. Finally each ICurve object consists of exactly two Point objects: the start and end point of this curve. In order to clarify the above description consider the nonconforming mesh in Fig. 4.2. Its associated mesh data structure is provided in table 4.2. Notice that the edges of a face are ordered in clock-wise direction. Similarly for curves, the index of the start point is always lower than the index of the endpoint. For fully discontinuous Galerkin settings the order doesn’t matter. However if the method is conforming in some parts of the mesh, the orientation of the edge is important to fix the orientation of the Nedelec

elements. Table 4.2 also shows that in 2D every face belongs to exactly one face segment

and that they both share the same edges.

Remark. Usually the IFace, IEdge and the Point objects have a direct analogue in

NGSolve’s MeshAccess interface. The corresponding ID can be obtained by calling GetNGID().

Remark. The edgesegments are constructed by calling

NonConformingMesh.ConstructEdgeSegments() in O(N log(N)) time where N 37

SLIDE 41

4. Aspects of the Implementation

Figure 4.1.: UML diagram of the mesh data structure. 38

SLIDE 42

4.2. Details of the imlementation. Figure 4.2.: Example of nonconforming mesh. Table 4.2 shows the associated mesh data structure. Faces FaceSegment Edges f0 fs0 (e0, e1, e2) f1 fs1 (e3, e4, e5) f2 fs2 (e6, e7, e8) Edge ID EdgeSegments Start-End Point e0 (es0) (p0, p1) e1 (es1) (p0, p3) e2 (es2, es3) (p1, p3) e3 (es4) (p1, p2) e4 (es2) (p1, p4) e5 (es5) (p2, p4) e6 (es5) (p2, p4) e7 (es3) (p3, p4) e8 (es6) (p2, p3) Table 4.2.: Mesh data structure for sample mesh in Fig. 4.2. 39

SLIDE 43

4. Aspects of the Implementation

is the number of edges. To achieve this the algorithm ”sorts” all curves according to their geometric parametrization (we use the term parametrization to refer to some parameters of the curve which describe it uniquely). After they have been sorted, all curves which could potentially intersect (have same parameterization) are neighbours. The edge segments are then constructed easily.

Remark. All the methods in the Surface<TRIG>, Curve<LINE> and Curve<CIRC_SEGM

have been marked as final (introduced in C++0x), i.e. all deriving classes cannot over- ride the implementation. This allows the compiler to resolve the virtual functions of these classes at compile time. Therefore efficient algorithms can be written by using static casts for specific Curve/Surface classes. On the other hand all the Curve resp. Surface objects have a common interface and can therefore easily be stored together in

ne container.
Remark. In principle the IFace, IFaceSegment, IEdge and IEdgeSegment classes

could also use composition instead of inheritance to provide the geometric methods from the Surface and Curve classes. E.g. the IEdge<LINE> class could internally store a reference to a Curve<LINE> object and delegate the respective methods to it. However this poses a rather big problem for the previously mentioned sorting algorithm. In this algorithm IEdge<LINE> objects are compared with each other. The comparison

perator is defined in the base class Curve<LINE> and needs to access the private

members of the other object. This is not possible when using composition.

Circle Segments

In principal the mesh data structure supports triangular IFace objects which have IEdge<CIRC_SEGM> objects as edges. Strictly speaking this doesn’t make sense because the IFace objects naturally define an element transform that maps from a reference triangle to another triangle. I.e. points which lie on the boundary of the reference triangle would not be mapped to points on the circle segment edges. Nevertheless we have used this approach to properly model domains with circle shaped

boundaries. Consequently the surface integrals are not analytically exact and contribute

an error of O(h2). Note that the flux integrals have a much smaller error: They are computed based on the FaceSegment<CIRC_SEGM> objects which in turn provide a element transformation which maps from a reference segment to the circle segment. We use this transformation together with an integration rule of order 10 which leads to an almost negligible error.

4.3. Variational formulation for moving systems

We have already discussed how the SWIP-formulation of the eddy current problem allows for an easy handling of non-matching grids. In this section we further explore the physical effects which appear if part of the system is in motion with respect to the other one and 40

SLIDE 44

4.3. Variational formulation for moving systems Figure 4.3.: The integration area over which the surface integrals are evaluated and the path along which the Flux integrals are calculated. how we can incorporate them into the variational formulation. We are going to base our discussion on the transformation properties which have been derived in Sec. 2.4. Let us first observe that the transformation matrix T which appears in relations (2.19) can be neglected for scalar problems. Indeed, a scalar quantity represents a vector with zero x − y coordinates. However the T-matrix acts only on these components. If we are dealing with a 2D vector problem we have to keep in mind that the implementation uses edge functions, i.e. T ˜ A is actually the quantity which is represented by basis functions and not ˜ A. Consequently, the magnetic field H doesn’t change in our implementation under translation and rotation, i.e. T ˜ H = H. Hence the problems (4.2) and (4.4) need no modification if part of the mesh is translated or rotated. In other words the convective terms which result from the Lorentz forces are automatically and correctly accounted for in the H − formulation. Consider now the TM-A potential formulation (4.3) with the primary, scalar unknown Ah,z. The correct transformation rule is given by lemma 2.3. However notice that the term grad

V · (0, 0, Ah,z)T

≡ 0 . Therefore ˜ Ah,z = Ah,z as well and the formulation (4.3) needs no modification. Finally we consider the TE-A potential formulation (4.1). Here the transformation law is given by 2.3 and reads as T ˜ A = A − T

t

TT grad (V · A). Using Eq. (A.11) we can simplify the last term as follows: grad (V · A) = −ω × A + (V · ∇)A + V × curl A. 41

SLIDE 45

4. Aspects of the Implementation

Next we assume that A ∈ R1(Th) and make use of Lemma 4.1. Let u ∈ R1(T), T ∈ Th and let V ∈ [L2(T)]3. Then there holds (V · ∇)u = −1 2V × curl u

Proof. We start from the representation formula u = λi · ∇λj − λj · ∇λi where λi, λj ∈

P1

3(T) (cf. [13, Eq. 5.47]). Straightforward manipulations then show that the assertion

holds. Therefore the transformation formula for A in (4.1) is T ˜ A = A + T

t

TT (ω × A) − T

t

TT (V × curl A). (4.5) Now consider the two general problems presented in Fig. 4.4. Here Ω2 is always at rest while Ω1 is either rotating or translating. We will assume that continuous basis functions are used in both subdomains such that the flux integrals vanish inside of Ω1 and Ω2. The flux is only non-zero at the boundaries resp. interfaces Γ1, Γ2 and Γ12.

(a) Translation problem. (b) Rotation problem

Figure 4.4.: Two possible eddy current problems that involve moving bodies. In order to calculate the fluxes over Γ12 we transform the vector potential A from the rest frame Ω2 to the moving frame Ω1. Notice that the tangential component of V × curl A vanishes at the interface Γ12 in both cases.

Variational formulation for translating bodies

We start from the variational formulation (4.1). Because ω = 0, we have for the tangential jump nF ×

T ˜

A1 − T ˜ A2

= nF ×
T ˜

A1 − A2

. This means that we don’t have to modify any terms on the left-hand side

42

SLIDE 46

4.3. Variational formulation for moving systems

f (4.1). On the right-hand side the source term stays exactly the same (T˜

ji = ji) and for the Neumann boundary term we have ˜ g = g. Indeed, ˜ g = 1

µ

˜ curl ˜ A = 1

µ curl A = g

because curl grad A = 0. Therefore the variational formulation (4.1) needs no modification as long as ω = 0.

Variational formulation for rotating bodies

Here the boundary Γ1 = ∅ and therefore the right-hand side of (4.1) stays exactly the same. However the left-hand side needs some

modification. First observe that the tangential jump across Γ12 is

nF ×

T ˜

A(i+1)

1

− T ˜ A(i+1)

2

= nF ×
T ˜

A(i+1)

1

− A(i+1)

2

− T

(i+1)·δt

TT (ω × A2(t))

,

= nF ×

T ˜

A(i+1)

1

− A(i+1)

2

− δt 2 ω × A(i+1) −δtT(i+1)

i

j=1

TT,(j)ω(j) × A(j)

2

− δt 2 T(i+1)(ω(0) × A(0)

2 ) + O(δt2)

  ,

where we have used the trapezoidal integration rule to approximate the integral over

time. Substituting this into the variational formulation (4.1) we get,

Find A(j)

h

∈ Vh := (Ah ∈ H(curl, Ω1) × H(curl, Ω2)) ∩ R1(Th), j = 1, . . . , N, such that for all A′

h ∈ Vh, all i = 0, . . . N − 1,

σA(i+1)

h

− A(i)

h

δt , A′

h

+ aSWIP

h,rot (A(i+1) h

, A′

h) =

ji,(i+1), A′

h

+
∂Ω

(n × A′

h) · g(i+1) + li rot(A′ h),

A(0)

h

= πh(A(t = 0)), where aSWIP

h,rot is a modified version of the bilinear form µaSWIP h,T

: aSWIP

h,rot (Ah, A′ h) =

Ω

1 µ curl2D,h Ah · curl2D,h A′

h

−

F∈Γ12
F

1

µ curl2D,h Ah

ω

·

A′

h

T − δt

2 nF ×

ω(i+1) × A′

h,2

−
F∈Γ12
F

1

µ curl2D,h A′

h

ω

·

[Ah]T − δt

2 nF ×

ω(i+1) × Ah,2
+
F∈Γ12

ηγµ,F hF

F
[Ah]T − nF ×

δt

2 ω(i+1) × Ah,2

·

A′

h

T − δt

2 nF ×

ω(i+1) × A′

h,2

.

(4.6) 43

SLIDE 47

4. Aspects of the Implementation

Note that this bilinear form is also symmetric. The linearform li

rot is given by

li

rot(A′ h) = −

F∈Γ12
F

1

µ curl2D,h A′

h

ω

·

 nF ×  δtT(i+1)

i

j=0

ǫjTT,(j) ω(j) × A(j)

h,2



  

+

F∈Γ12

ηγµ,F hF

F

 nF ×  δtT(i+1)

i

j=0

ǫjTT,(j) ω(j) × A(j)

h,2



   · A′

h

T − δt

2 nF ×

ω(i+1) × A′

h,2

,

where ǫ0 = 1/2 and ǫj>0 = 1.

Remark. We have not only transformed the trial functions Ah but also the test functions

A′

h! In theory this is not needed but it gives a symmetric matrix which forms, together

with the mass matrix, a symmetric-positive-definite matrix. In practice it was also

bserved that the transmission conditions are much better fulfilled if the test functions

are transformed in the above way.

Remark. We have used the trapezoidal rule to approximate the time integral which

appears in the transformation rule (4.5). In the actual implementation the values of the integral are interpolated to the midpoints of every edge and are integrated at this location.

Remark. We have seen that the eddy current problem in the H-formulation and the time

gauged potential formulation need no modification for translating bodies (This can easily be extended to 3D). In other words just moving the mesh already takes into account the convective terms. A similar statement was proven by Kurz et al. [10, Section II] for the Coulomb gauged eddy current formulation. 44

SLIDE 48

5. Results

After having discussed the theoretical aspects of the DG formulation and after we have described the implementation on the computer we finally present some numerical results. In particular we will discuss three examples that shed light on different aspects of the formulations. In Sec. 5.1 we will consider a radial problem for which an analytical solution is known. We can therefore investigate the theoretical convergence bounds which have been proven in chapter 3. Unfortunately this problem doesn’t contain any convective terms by con-

struction. Constructing an analytical solution in which the convective terms play a role

is very hard. To overcome this problem we resort to comparing the TE-time gauged formulation (4.1) with the TE-H formulation. They both describe the same physical problem so similar results are expected. Section 5.2 studies this for translational motion whereas Sec. 5.3 does the same for a rotating body.

5.1. Convergence to analytical solution

We consider the following problem in polar coordinates on the domain Ω = {(x, y) ∈ R2

x2+

y2 < 1}: ∂Hz ∂t − div grad Hz = J0(λ1r)

λ2

1

2 t2 + t

,

(5.1a) Hz(r, t = 0) = 0, (5.1b) Hz(r = 1, t) = 0, (5.1c) where J0 is the first Bessel function of the first kind and λ1 is it’s first root (λ1 ≈ 2.4048255). One can show that this IBVP posses the unique solution Hz(r, t) = 1 2t2J0(λ1r).

SLIDE 49

5. Results

The problem (5.1) can be discretized straightforwardly and gives the TE-H formulation 4.2 introduced in chapter 4 (µ = σ = 1). We split the domain Ω into two subdomains Ω1 = {(x, y) ∈ R2

x2 + y2 < 1/4}, Ω2 = Ω\Ω1. The inner circle is then additionally

being rotated (circle segments have been introduced along Ω1 ∩ Ω2) with ω = 20 rad/s. Figure 5.1a shows the two subdomains whereas Fig. 5.1b presents the calculated solution

f the problem at t = 0.05 s .

(a) The subdomains Ω1 and Ω2 (b) The solution at time t = 0.05 sec(δt = 0.0005, h = 0.0361371)

Figure 5.1.: Discretization of the IBVP 5.1

Convergence to the exact solution

From Theorem 3.21 in chapter 3 we expect to get O(h) respectively O(δt) convergence to the analytical solution. In order to prove this we have carried out a number of numerical experiments with different mesh sizes and

timesteps. In particular we study the error in the L2 norm and the (modified) H curl

norm A − Ahcurl :=

A(N) − A(N)

h

L2(Ω) +
Cstabδt

N

i=1
A(i) − A(i+1)

h

2

SWIP

1/2

. The error in both norms has been approximated by an integration rule of order 5 on each triangle and of order 10 on each edge. The total simulation time was fixed to 0.05 s and only the timestep δt and the mesh size h were varied.

Figs. 5.2a and 5.2b show that we get convergence in both norms if δt as well as h are

scaled properly. This is exactly what we expect from Thm. 3.21. Next take a look at

Figs. 5.3a and 5.3b which show the error in both norms for fixed δt respectively fixed h.

We see that the measured rate of convergence in terms of δt is slightly better than the theoretical prediction whereas the convergence in h is slightly worse for the H curl-norm. 46

SLIDE 50

5.1. Convergence to analytical solution

(a) L2 error (b) H curl error

Figure 5.2.: Error vs. δt and h.

(a) δt convergence, h = 0.0180874 (b) h convergence, δt = 2.5e − 4

Figure 5.3.: Convergence plots for h and δt 47

SLIDE 51

5. Results
Remark. Di Pietro and Ern [6, Sec. 4.7.5] show that the theoretical rate of convergence

is O(h2) in the L2 norm. This estimate is optimal but is not achieved in this example because of the approximation error on the circle segments (see 4.2.1).

5.2. Translational motion

The two previous test cases have shown that the TE-H formulation (4.2) produces the expected results in the presence of sliding interfaces. In this section we will compare this formulation with the TE-A formulation and measure the (relative) error between the two approximate solutions. Figure 5.4.: Dimensions of the translational test case. We consider the setting depicted in figure 5.4 The bottom rectangular body is at rest while the top one slides along its edge with velocity V = 1 m/s. The material coefficients are all set to unity (µ = σ = 1) and we use homogeneous Neumann boundary conditions for the TE-A formulation. Consequently we apply homogeneous Dirichlet boundary conditions for the H-field. The system is excited by the impressed source current ji =

      

(0, 0)T for x ≤ 1 or x ≥ 3 (0, x − 1)T for 1 < x < 2 (0, 3 − x)T for 2 ≤ x < 3 . The total simulation time is tEnd = 1 s and we start from the initial condition A = Hz = 0 Am−1. Figure 5.5 shows the solutions shortly after two timesteps and at the at the final time. We can clearly see that the motion of the top body influences the shape of the H-field and that both formulations yield similar results. We now turn to a more quantitative analysis and study the error in terms of meshwidth and timestep size. The error Hz − curl2D,h AL2(Ω) is measured slightly different for the two cases: h-convergence We simulate both systems with the same timestep δt and the same mesh size and measure the L2 error between them for successive refinements of h (both formulations are refined). 48

SLIDE 52

5.2. Translational motion

(a) A-solution (curl2D,h A shown in the background) (b) H-solution

Figure 5.5.: Solutions of the translational test case at t = 0.01 s (top) respectively t = 1 s (bottom). δt = 5 · 10−3 s, h = 0.1346 m. 49

SLIDE 53

5. Results

δt convergence We compute a (reference) solution of the H-formulation with timestep δt = 10−3 s and meshwidth h = 0.0336515 m. The temporal gauged potential formulation is simulated at the same mesh width but we vary the size of the

timestep. The L2 error is then calculated with respect to the H-formulation.
Remark. In the limit of very small h it doesn’t matter how we measure δt convergence:

Taking one model as a reference and while refining the other produces the same rate of convergence as refining both models together: For both approaches we essentially plot log(δt + ε(h)) + const = log(δt) + ε(h)/δt + O(ε(h)2) + const where ε(h) → 0 as h → 0. However if we take one solution as a reference ε(h) is typically smaller than if we refine both systems together. Figure 5.6 shows the convergence in terms of h and δt. The measured rates are 0.79 for h and 0.922 for δt. This is really close to the values of Thm. 3.21.

(a) h convergence, δt = 5 · 10−3. (b) dt convergence, h = 0.03365

Figure 5.6.: Convergence for the translational test case.

5.3. Rotational motion

Finally we want to compare the H-model with the temporal gauged potential formulation for a rotating body. In contrast to the previous simulation additional transmission terms appear in the variational formulation of A (see Eq. (4.6)). We consider the setting depicted in Fig. 5.7: The inner circle is rotating at a speed of ω = π rad/s and the total simulation time is 1 s. The material coefficients are µ = σ = 1 in both sub-domains. We consider the TE formulations (4.2) and (4.1) for the H-field respectively the A-

field. We use homogeneous Neumann boundary conditions for A and correspondingly

50

SLIDE 54

5.3. Rotational motion Figure 5.7.: The setup of the numerical experiment considered in this section. homogeneous Dirichlet boundary conditions for H. On the right hand side we impose the source current ji = (4y, −2x). Both simulations are carried out at the same time and we start from the initial condition A = Hz = 0. Figure 5.8 shows the evolution of the A field (first column) and the Hz field (last column) over time. In addition the solution which is obtained for the ”simple” tangential boundary conditions is also included in the middle column (i.e. we used aSWIP

h,T

from Sec. 4.1 instead of aSWIP

h,rot ) . At the beginning the system equilibrates (0 ≤ t 0.08). In this

phase the curlh A respectively Hz of all three solutions look qualitatively similar but as time evolves only the H field and the curlh A (with additional transmission conditions) look alike. The middle column converges clearly to the wrong solution. It’s interesting to see that the curlh A field (with additional transmission condition) is somewhat time- periodic; after one full rotation (t = 1 s) the vector field looks very similar to the vector field at t = 0.05 s. Indeed we observe that the A vector field periodically rotates up to the position at t ≈ 0.95 s and then switches back to the initial position. The period is approximately 1 second. Figure 5.9 shows the convergence rates as δt resp. h decrease. We measure the error in the same way as for the previous test case: For h-convergence we refine both systems at the same time. In order to measure the δt convergence we use the H-formulation as a reference solution (δtH = 10−3) and measure the error to the A-formulation for different sizes of the timestep δtA. We see clearly that we get convergence in both variables. Although this is only the relative error between the H-formulation and the A-formulation we can still conclude that they approximate the analytical solution asymptotically at a rate of approximately 0.5 resp. 0.7 (assuming that neither of the two converges to the wrong solution). It remains unclear as to why the observed rates are not of order 1. If the discrete compactness property (which was used in the derivation of Thm. 3.21) can be extended to the more general jump-seminorm appearing in the aSWIP

h,rot

(Eq. (4.6)) it should be 51

SLIDE 55

5. Results

t = 0.05 t = 0.5 t = 0.95 t = 1.04

(a) A (vectors), curlh A (background) if transmission condition (4.5) is used (b) A (vectors), curlh A (background) if only tangential continuity is enforced. (c) Hz field (background)

Figure 5.8.: Different solutions of the rotational test case. (h = 0.101761) 52

SLIDE 56

5.3. Rotational motion

(a) δt convergence for different mesh sizes (b) h convergence for δt = 0.001 sec.

Figure 5.9.: Convergence of the relative error Hz − curl2D,h AL2(Ω) at tend = 1 (one full rotation) possible to proof a result similar to Thm. 3.21. In fact the author has considered a much simpler test case with a constant, impressed current ji on the right-hand side. With this setup linear O(δt) convergence was observed. 53

SLIDE 57

6. Conclusion and Outlook

Maxwell’s equations are invariant under the Lorentz transformation if the moving frame

f reference is travelling at uniform speed. This is not true for rotation because already

small accelerations have a significant effect (see [19], [10]). However, in the quasi-static limit of slowly varying electric fields we were able to prove the invariance of the Eddy Current model under solid body motion (x = T(t)˜ x + r(t)). This invariance justifies a Lagrangian description. In order to deal with the non-matching interfaces we have adopted the Discontinuous Galerkin technique and were able to give a proof for the simple case of matching grids and (sufficiently smooth) H(curl, Ω) conforming solutions. This result was verified with the 2D, scalar H-formulation for non-matching grids which suggests that the convergence theorem also holds on nonconforming meshes. In fact the theorem could be easily extended if the discrete compactness property ([12]) could be proven for non-matching meshes. Furthermore it would make sense to test the convergence theorem in a 3D setting. Because the transformation law T ˜ B = B holds for the magnetic induction the variational formulation needs no change when treating problems with moving bodies (if edge functions are used!). However a special transmission condition must be used for the A temporal gauged potential formulation in rotational settings; the transmission conditions from the steady case do not suffice! This was verified with numerical experiments for the TE-formulations: For the translational test case we have measured the expected rates of convergence. For the rotational setting we would expect a rate of O(δt) but we have measured the convergence rate O(δt

1/2). It is not clear what causes the difference

and further investigations are needed. For the numerical experiments we have developed a (rather complex) mesh data structure to represent the non-conforming nature of the mesh. The algorithm to detect edge segments worked without problems for very fine meshes in O(N log(N)) time. It is unclear if the algorithm’s idea can be extended in an efficient way to 3D, respectively 2D Face segments.

SLIDE 58

A. Maxwell’s equations in a moving frame of

reference

Consider a body which, from the view of the laboratory system, is moving in a direction r(t) rotating around an axis ω. Then we can establish a body-centered coordinate system ˜ x such that the laboratory coordinates x are given by, x = T(t)˜ x + r(t), (A.1a) t = ˜ t (A.1b) where T is a (orthogonal) rotation matrix only depending on time. The velocity of a fixed point ˜ x in the moving system is given by (in laboratory coordinates): V := ∂x ∂t = ∂T ∂t ˜ x + ∂r(t) ∂t , = ∂T ∂t TT (x − r(t)) + ∂r(t) ∂t . On the other hand we can express the velocity also in terms of the angular velocity ω: V = ω × (x − r(t)) + ∂r(t) ∂t Comparing the two expressions we see that for any vector a ∈ R3: ∂T ∂t TT a = Wa = ω × a with W =

  

−ω3 ω2 ω3 −ω1 −ω2 ω1

   .

(A.2) Our goal is to state the reduced Maxwell equations in in terms of ˜ x and it’s derivatives. In the laboratory system the equations are, div B = 0 (A.3a) curl E + ∂B ∂t = 0 (A.3b) curl H = ji + jf. (A.3c)

SLIDE 59

A. Maxwell’s equations in a moving frame of reference
Remark. We have left out the electric Gauss law on purpose because there seems to be

no meaningful transformation of the electric induction D under transformation (A.1a). This is not a problem because this equation is completely decoupled from the rest and is not needed to formulate the eddy current problem. Theorem A.1. Assuming the constitutive relations of linear materials,the system of equations (A.3) is invariant under the transformation (A.1a) if the field quantities transform in the following form (cf. [10]): T˜ E = E + V × B T ˜ B = B T˜ ji = ji T ˜ H = H T˜ jf = jf.

Proof. First notice that the del-operator in the body’s coordinate system can be written

as: ˜ ∇ = Dx˜ x ∇ = TT ∇. Therefore Eq. (A.3a) in the moving frame of reference can be rewritten as ˜ ∇ · ˜ B =

TT ∇

T

TT B

= ∇T TTT B = ∇ · B = 0.

This means that ˜ div ˜ B = 0 holds and therefore the magnetic Gauss law holds also in the moving system. In order to prove the same for Faraday’s law of induction (A.3b) we first notice that because (Ma) × (Mb) = (det M)M−T (a × b), we have ˜ ∇ × A = (det TT )T−1 (∇ × (TA)) = TT (∇ × (TA)) . (A.4) Similarly we can express ∂ ˜

B ∂˜ t as

∂ ˜ B

x(˜

x, ˜ t), t(˜ t)

∂˜

t = ∂TT B ∂˜ t = ∂TT ∂t B + TT ∂B ∂˜ t

(A.2)

= (WT)T B + TT ∂B

x(˜

x, ˜ t), t(˜ t)

∂˜

t = −TT (ω × B) + TT

∂B

∂t +

∂x

∂˜ t · ∇

B
.

Notice that ∂x

∂˜ t = V and use the vector calculus identity ∇ × (V × B) = V(∇ · B) − B(∇ ·

V) + (B · ∇)V − (V · ∇)B. The first term is obviously zero because magnetic charges do not exist, the second term simplifies to −B(∇ · V) = −B

∇ ·
ω × (x − r(t)) + ∂r

∂t

=

−B (∇ · (ω × x)) = −x · (∇ × ω) + ω · (∇ × x) = 0 and for the third term we get: (B · ∇)V = (B · ∇)(ω × x) = (B · ∇)(Wx) = WB = ω × B. (A.5) 56

SLIDE 60

Therefore, ∂ ˜ B

x(˜

x, ˜ t), t(˜ t)

∂˜

t = TT

−ω × B + ∂B

∂t + ∇ × (V × B) + ω × B

= TT

∂B

∂t − ∇ × (V × B)

.

(A.6) Now substituting Eqs. (A.4), (A.6) into Faraday’s law and using the transformation rules of Thm. A.1 we obtain ˜ curl˜ E + ∂ ˜ B ∂˜ t = TT

curl E + curl (V × B) + ∂B

∂t − curl (V × B)

= TT
curl E + ∂B

∂t

= 0.

I.e. Faraday’s law of induction is also invariant under the transformation (A.1a). For Ampères law (A.3c) we use Eq. (A.4) to get, ˜ curl ˜ H − ˜ ji − ˜ jf = TT curl H − ji − jf = 0.

Remark. It is well known that the full Maxwell equations are invariant under the Lorentz

transformation (see e.g. [16, chp. 22]). In the quasistatic limit (i.e. c → ∞) the lorentz transformation reduces to a Galilean transform and the fields are transformed as described by Thm. A.1. Thm. A.1 is however more general than this as it considers general solid body motion which is generally not a Galilean transform. (T(t) as well as r(t) can depend on time!) This is however only possible in the quasistatic limit. For high frequencies where the displacement current can generally not be neglected, the acceleration of a rotating frame

f reference can have significant effects on the field quantities.

Generalized Ohm’s law

It should be stressed that Ohm’s law, jf = σE is only valid if the underlying material is at rest. Let us assume that a body is moving according to

Eq. (A.1a). Then obviously ˜

jf = σ ˜

E. Therefore, by virtue of Thm. A.1, the relation

jf = σ(E + V × B) (A.7) holds in the laboratory’s system. The second term on the right-hand side accommodates for the Lorentz forces acting on the moving charges. This is consistent with a "physical" derivation of Ohm’s law for a rotating body (cf. [11]). 57

SLIDE 61

A. Maxwell’s equations in a moving frame of reference

Transformation rules for potential formulation

In the potential formulation of the eddy current problem one uses a a vector and scalar potential such that, B = curl A∗ E = − grad ϕ∗ − ∂A∗ ∂t . (A.8) The question arises now as to how these potentials transform under Eq. (A.1a): Corollary A.2. Under the assumptions of Theorem A.1, the equations (A.8) are invariant under transformation (A.1a) if the field quantities transform as follows: T˜ A∗ = A∗ ˜ ϕ∗ = ϕ∗ − V · A∗ Together with Thm. A.1 this implies that if A∗ and ϕ∗ are solutions to Eqs. (A.3) then ˜ A∗ and ˜ ϕ∗ are the respective solutions in the moving frame.

Proof. By using formula (A.4), it is easily seen that ˜

B = ˜ curl ˜ A∗ = TT curl A∗ = TT B. This is exactly the transformation law of B and therefore the first equation of (A.8) is invariant if T ˜ A∗ = A. In order to prove that E = − grad ϕ∗ − ∂A∗

∂t

is invariant, we start from the laboratory frame and transform the equation to express it in terms of the body’s coordinates: E = − grad ϕ∗ − ∂A∗ ∂t ⇔ T˜ E − V × curl A∗ = − grad ˜ ϕ∗ − grad (V · A∗) − ∂T ˜ A∗ ∂t (A.9) Now use that, grad (V · A∗) = (A∗ · ∇)V + (V · ∇)A∗ + A∗ × (∇ × V) + V × (∇ × A∗)

(A.5)

= ω × A∗ + (V · ∇)A∗ + A∗ × (∇ × (ω × x)) + V × (∇ × A∗) ⇒∇ × (ω × x) = ω(div x) − x(∇ · ω) + (x · ∇)ω − (ω · ∇)x = 3ω − ω = 2ω = −ω × A∗ + (V · ∇)A∗ + V × (∇ × A∗), (A.10) and, ∂T ˜ A∗ ∂˜ t = ∂T ˜ A∗ ∂t + (V · ∇)(T ˜ A∗) ⇔ ∂T ∂t ˜ A∗ + T∂ ˜ A∗ ∂t = ∂T ˜ A∗ ∂t + (V · ∇)(T ˜ A∗) ⇔ ∂T ˜ A∗ ∂t = ω × A∗ + T∂ ˜ A∗ ∂˜ t − (V · ∇)A∗. (A.11) 58

SLIDE 62

Therefore Eq. (A.9) can be rewritten as, T˜ E − V × curl A∗ = −T ˜ grad ˜ ϕ∗ + ω × A∗ + (V · ∇)A∗ − V × curl A∗ − ω × A∗− T∂ ˜ A∗ ∂˜ t − (V · ∇)A∗ ⇔ T˜ E = T

−

˜ grad ˜ ϕ∗ − ∂ ˜ A∗ ∂˜ t

59

SLIDE 63

B. Documentation of the code

This chapter together with together with the Doxygen documentation provides a detailed description of the code which has been written as a part of this thesis. The source folder of this thesis contains five folders:

NonConforming In this folder lies the NonConforming library. It provides almost all

the reusable functionality that is needed to solve the eddy current problem on non-conforming meshes. The Doxygen documentation of this library is containted in the file refman.pdf which lies in the NonConforming folder. It gives a very detailed description of the actual classes.

NumProcs In this folder all the numerical experiments which have been presented in

chapter 5 lie in this folder. It is essentially a collection of custom NumProc classes which can be loaded at runtime by NGSolve/Netgen.

NonConformingTest All the numerical tests which have been written with the unit-

testing framework CUTE [17] lie in this folder.

Latex contains the sources of this document. NGSolve Contains the sources of NGSolve 4.9. The author designed the NonConforming

library such that NGSolve’s and Netgens code needed as little modification as

possible. However a few minor changes were necessary to make everything working

and they are included in this folder.

Remark. Everything was built on NGSolve 4.9 and Netgen 4.9.14-RC.

SLIDE 64

C. Acknowledgements

I would like to thank Prof. Hiptmair for the interesting discussions and the many helpful comments which lay the foundation for this thesis.

SLIDE 65

Bibliography

[1] A high order discontinuous galerkin finite element solver for the incompressible navier–stokes equations. Computers & Fluids, 46(1):224 – 230, 2011. (Cited on page 2.) [2] R Albanese and G Rubinacci. Formulation of the eddy-current problem. IEE Proceedings A (Physical Science, Measurement and Instrumentation, Management and Education), 137(1):16–22, 1990. (Cited on pages 3 and 32.) [3] Alain Bossavit. Whitney forms: a class of finite elements for three-dimensional computations in electromagnetism. Physical Science, Measurement and Instrumen- tation, Management and Education-Reviews, IEE Proceedings A, 135(8):493–500,

1988. (Cited on pages 11 and 27.)

[4] Annalisa Buffa, Yvon Maday, and Francesca Rapetti. A sliding mesh-mortar method for a two dimensional eddy currents model of electric engines. ESAIM: Mathematical Modelling and Numerical Analysis, 35(02):191–228, 2001. (Cited on pages 2 and 32.) [5] Emmanuel Creusé and Serge Nicaise. Discrete compactness for a discontinuous galerkin approximation of maxwell’s system. ESAIM: Mathematical Modelling and Numerical Analysis, 40(02):413–430, 2006. (Cited on pages 17 and 20.) [6] Daniele Antonio Di Pietro and Alexandre Ern. Mathematical aspects of discontinuous Galerkin methods, volume 69. Springer, 2011. (Cited on pages 11, 12, 13, 14, 15, 16, 17, 18, 20, 30, 31, 34 and 48.) [7] Esteban Ferrer. A high order Discontinuous Galerkin - Fourier incompressible 3D Navier-Stokes solver with rotating sliding meshes for simulating cross-flow turbines. PhD thesis, University of Oxford, 2012. (Cited on page 2.) [8] B. Flemisch, Y. Maday, F. Rapetti, and B.I. Wohlmuth. Coupling scalar and vector potentials on nonmatching grids for eddy currents in a moving conductor. Journal

f Computational and Applied Mathematics, 168(1–2):191 – 205, 2004. ISSN 0377-

0427. doi: 10.1016/j.cam.2003.05.017. URL http://www.sciencedirect. com/science/article/pii/S0377042703009774. (Cited on page 2.) [9] John David Jackson. Classical Electrodynamics. Wiley, third edition, 1998. ISBN

SLIDE 66

BIBLIOGRAPHY 978-0471309321. (Cited on page 32.) [10] S. Kurz, J. Fetzer, G. Lehner, and W.M. Rucker. A novel formulation for 3d eddy current problems with moving bodies using a lagrangian description and bem-fem coupling. Magnetics, IEEE Transactions on, 34(5):3068–3073, 1998. (Cited on pages 44, 54 and 56.) [11] L.D. Landau, E.M. Lifsic, J.B. Sykes, J.S. Bell, MJ Kearsley, and L.P. Pitaevskii. Electrodynamics of continuous media, volume 364. Pergamon press Oxford, 1960. (Cited on page 57.) [12] P Monk and L Demkowicz. Discrete compactness and the approximation of maxwell’s equations in R3. Math. Comp, 70(234):507–523, 2001. (Cited on pages 13, 17, 31 and 54.) [13] Peter Monk. Finite Element Methods for Maxwell’s Equations, Numerical Math- ematics and Scientific Computation. Oxford University Press, New York, 2003. (Cited on pages 11, 17, 20, 21, 27, 31, 34 and 42.) [14] F Rapetti, E Bouillault, L Santandrea, A Buffa, Y Maday, and A Razek. Calculation

f eddy currents with edge elements on non-matching grids in moving structures.

Magnetics, IEEE Transactions on, 36(4):1351–1355, 2000. (Cited on page 2.) [15] Francesca Rapetti, Yvon Maday, Frédéric Bouillault, and Adel Razek. Eddy-current calculations in three-dimensional moving structures. Magnetics, IEEE Transactions

n, 38(2):613–616, 2002. (Cited on page 2.)

[16] J.R. Reitz, F.J. Milford, and R.W. Christy. Foundations of electromagnetic theory. Fourth edition, 1995. (Cited on page 57.) [17] Peter Sommerlad and Emanuel Graf. Cute: C++ unit testing easier. In Companion to the 22nd ACM SIGPLAN conference on Object-oriented programming systems and applications companion, pages 783–784. ACM, 2007. (Cited on pages 35 and 60.) [18] J. Van Bladel. A discussion of helmholtz’theorem. Electromagnetics, 13(1):95–110,

1993. (Cited on page 4.)

[19] DavidL. Webster and RobertC. Whitten. Which electromagnetic equations apply in rotating coordinates? Astrophysics and Space Science, 24(2):323–333, 1973. ISSN 0004-640X. URL http://dx.doi.org/10.1007/BF02637159. (Cited

n page 54.)

[20] Sabine Zaglmayr. High Order Finite Elements for Electromagnetic Field Computa-

tion. PhD thesis, Johannes Kepler Universität Linz, 2006. (Cited on page 5.)

63

Sliding Interfaces for Eddy Current Simulations Raffael Casagrande - - PDF document

Sliding Interfaces for Eddy Current Simulations

Raffael Casagrande

Master Thesis Supervisor:

Contents

2.1. Coulomb gauged formulation

2.2. Temporal gauged formulation

2.3. H-formulation

2.4. The eddy current problem in a moving, solid body

3.1. Preliminaries

3.2. DG formulation of Eddy Current problem

4.1. Reduction to 2D

4.2. Details of the imlementation.

4.3. Variational formulation for moving systems

5.1. Convergence to analytical solution

5.2. Translational motion

5.3. Rotational motion

reference

Bibliography