FEM Implementierung von Spulen Modellen f ur Transformator ... · Transformator Simulationen FEM...

Diese Arbeit wurde vorgelegt amLehrstuhl fur Mathematik (MathCCES)

FEM Implementierung von Spulen Modellen furTransformator Simulationen

FEM Implementation of Winding Models forTransformer Simulations

BachelorarbeitComputational Engineering Science

Marz 2018

Vorgelegt von Christian GellerPresented by Mauerstrasse 40, 52064 Aachen

Matrikelnummer: [email protected]

Erstprufer Prof. Dr. Manuel TorrilhonFirst examiner Lehrstuhl fur Mathematik (MathCCES)

RWTH Aachen University

Zweitprufer Prof. Dr. Benjamin StammSecond examiner Lehrstuhl fur Mathematik (MathCCES)

RWTH Aachen University

Koreferent Dr. Jorg OstrowskiCo-supervisor ABB Corporate Research, Power Devices Simulations

ABB Schweiz AG

Eigenstandigkeitserklarung

Hiermit versichere ich, dass ich diese Bachelorarbeit selbstandig verfasst und keineanderen als die angegebenen Quellen und Hilfsmittel benutzt habe. Die Stellen meinerArbeit, die dem Wortlaut oder dem Sinn nach anderen Werken entnommen sind, habeich in jedem Fall unter Angabe der Quelle als Entlehnung kenntlich gemacht. Dasselbegilt sinngemaß fur Tabellen und Abbildungen. Diese Arbeit hat in dieser oder einerahnlichen Form noch nicht im Rahmen einer anderen Prufung vorgelegen.

Aachen, 29. Marz 2018

Christian Geller

Contents

List of Figures IV

List of Tables V

1 Introduction 1

2 Robust Full Maxwell 52.1 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Numerical Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . 82.3 Solving Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

3 Excitation of Stranded Windings 133.1 Current Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.1.1 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . 143.1.2 Numerical Discretization . . . . . . . . . . . . . . . . . . . . . . 153.1.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.2 Voltage Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173.2.1 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . 173.2.2 Numerical Discretization . . . . . . . . . . . . . . . . . . . . . . 193.2.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.3 Calculation of Winding Tangential eϕ . . . . . . . . . . . . . . . . . . . 203.3.1 Failure of Analytical eϕ . . . . . . . . . . . . . . . . . . . . . . 213.3.2 Discrete divergence-free eϕ . . . . . . . . . . . . . . . . . . . . . 22

3.4 Postprocessing Calculations . . . . . . . . . . . . . . . . . . . . . . . . 26

4 Excitation of Massive Windings 294.1 Current Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 294.2 Voltage Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324.3 Postprocesssing Calculations . . . . . . . . . . . . . . . . . . . . . . . . 32

5 Validation 345.1 Verification of the Winding Tangential Calculation . . . . . . . . . . . 34

5.1.1 Analytical Solution . . . . . . . . . . . . . . . . . . . . . . . . . 345.1.2 Compare Analytical and Numerical Solution of ϕ0 and eϕ . . . 35

5.2 Internal Comparison With Existing Solver . . . . . . . . . . . . . . . . 375.2.1 Winding Model With Different Simulation Settings . . . . . . . 375.2.2 Compare Circular and Spiral Winding Model . . . . . . . . . . 39

5.3 External Comparison With Commercial Solver . . . . . . . . . . . . . . 405.3.1 Massive Winding . . . . . . . . . . . . . . . . . . . . . . . . . . 405.3.2 Realistic Transformer . . . . . . . . . . . . . . . . . . . . . . . . 42

6 Conclusion 46

Glossary 48

References 51

List of Figures

1 Spiral, Circular and Cylindric Model for six Turns . . . . . . . . . . . . 22 Spiral and Cylindric Model for a Three Layered Coil With 3000 Turns . 33 3D Tetrahedral Element . . . . . . . . . . . . . . . . . . . . . . . . . . 94 Code for Right Hand Side Calculation . . . . . . . . . . . . . . . . . . 175 Code for Quadrature

∫Teti

eϕ ·A dx . . . . . . . . . . . . . . . . . . . 206 Abstract Cross-Section of a Winding . . . . . . . . . . . . . . . . . . . 237 2D Discretized Cross-Section . . . . . . . . . . . . . . . . . . . . . . . . 258 Flow Chart of Implemented Structure in HADAPT . . . . . . . . . . . 279 Potential ϕ for Massive Voltage Excitation . . . . . . . . . . . . . . . . 3010 Results of the winding potential ϕ0 . . . . . . . . . . . . . . . . . . . . 3511 Results of the winding tangential eϕ . . . . . . . . . . . . . . . . . . . . 3612 Relative Error in L2 Norm in Log-log Plot of eϕ . . . . . . . . . . . . . 3613 Current Density j of Internal Comparison . . . . . . . . . . . . . . . . . 3814 Detailed Current Density j . . . . . . . . . . . . . . . . . . . . . . . . . 3815 Geometry of Circular and Spiral Model . . . . . . . . . . . . . . . . . . 3916 Current Density j on Wall for 50 Hz . . . . . . . . . . . . . . . . . . . . 4017 Current Density Re(j)z on Cross-Section for Different Frequencies . . . 4118 Current Density on Cross-Section for 500 Hz and Both Solvers . . . . . 4219 Simulation Setup for Transformer Model . . . . . . . . . . . . . . . . . 4320 Current Density Re(j)x on Windings for HADAPT and MagNet . . . . 4321 Current Density j on Housing for HADAPT and MagNet . . . . . . . . 4422 Magnetic Field B on Core for HADAPT and MagNet . . . . . . . . . . 45

List of Tables

1 Simulation Setup for Internal Comparison . . . . . . . . . . . . . . . . 372 Induced Losses in Transformer Model . . . . . . . . . . . . . . . . . . . 45

1 Introduction

1 Introduction

In general, the simulation of power devices and particularly of multi-phase transform-ers is essential in research and industrial companies. More precisely electromagneticproperties, like the electric or magnetic field, and further ohmic losses or inducedforces are important for the design of electronic devices. Hence these properties needto be simulated in the entire frequency domain. In the case of transformer simulationscoils with thousands of turns excite the system with a prescribed coil excitation. Theoptimization and simulation of, e.g., these transformer coils is the aim of this thesis.

State of the ArtAll electromagnetic properties can be referred to the general, analytical Maxwell equa-tions (1), so that the solving of these equations is required for the simulation of electricphenomena. A numerical formulation is derivated by [20] based on the general finiteelement method, which is observed in [4]. The more specific numerical approachesof solving the Maxwell equations are covered by [1] or [6]. Sometimes also the timedependency of the electromagnetic effects is considered in transient numerical compu-tations as in [18], but in the general, all properties are considered time harmonic aftera short calibration time as described in [14].

[7] derives a particular formulation of the Maxwell equations for the entire frequencydomain. Especially for the low-frequency spectrum, a robust formulation of the equa-tions is required for a solvable equation system, which leads to two introduced complexpotentials ϕ and A forming the solution variables. The mathematical equation systemcan be discretized with second order, curved tetrahedral and further solved with thenumerical finite element method. [7]

This robust, mathematical formulation of the electromagnetic problem is embeddedinside of a FEM solver, called HADAPT. For a given mesh discretization and appliedboundary conditions for voltages or currents the existing solver computes the robustformulated equation system and additionally all electromagnetic properties [8].

MotivationIn electronic power devices and especially in the case of multiphase transformers mul-tiple coils with a huge amount of turns exist, where a typical three-phase industrialtransformer contains six coils with several thousand turns each. The simulation ofthose coils is significant and already possible with the existing solver HADAPT. Hencea precisely meshed coil with spiral turns can be used in a simulation with an appliedcurrent or voltage as a boundary condition. As an example, a coil with six spiral turnsis depicted in Figure 1a.

In a first step, the spiral coil can be simplified by using circular instead of helical

1

1 Introduction

turns. This simplification leads to circular closed turns and allows the definition of theprescribed value separated for each turn inside of the domain. The simplification ofthe spiral coil is depicted in Figure 1b.

Because of the thin structure and vast amounts of turns, the meshing process of coilscan lead to problems or even be not feasible. Besides, the solving procedure itself isexpensive for a finely-meshed geometry so that a further simplification is needed forthose circular coils.

(a) spiral model(b) circular model

(c) cylinder model

Figure 1: Spiral, Circular and Cylindric Model for six Turns

In general, the structure and form of every spiral turn and therefore also the distribu-tion of turns is uniform along a coil, which leads to an equal current density in everysingle turn. In a simplified model, the coil can be pooled to one single hollow cylinder.The prescribed current or voltage is summed up overall turns and applied to the entirehollow cylinder. In addition, the meshing process becomes more viable because onlyone component without small edges or surfaces has to be meshed. The dimensions of

2

1 Introduction

this hollow cylinder are equal to the original coil because the new simplified cylinderrepresents the same amount of turns. The respective, simplified, cylindrical modelfrom the example above is visualized in Figure 1c.

A more complex model is a three-phase layered coil with 1000 turns per layer. For asimplified mesh generation and solubility, this model can also be expressed with the newapproach. Figure 2a and 2b show the physical spiral coil and the simplification withthree hollow cylinders for the different phase layers. Besides a less expensive meshing,the usual connections to the airbox for defining the excitations can be neglected.

(a) spiral model (b) cylinder model

Figure 2: Spiral and Cylindric Model for a Three Layered Coil With 3000 Turns

The new winding model, a simplification with circular hollow cylinders, is establishedand implemented for different winding types in the following thesis.

Types of WindingsThe first distinction of windings is the type of excitation. A primary coil always hasa prescribed excitation, which can be either a voltage or a current. Because of therelation between voltage and current, both of the two excitation types can be appliedand lead to the same physical results. For cylinders, which represent more than oneturn, the total applied current is the sum over all currents in the turns, which isIcylinder = Icoil ·Nturns for a uniform turn set up. In comparison the prescribed voltagefor several turns end up with Ucylinder = Ucoil/Nturns, due to the shortened perimeter ofthe voltage in a cylinder.

A prescribed excitation leads due to induction phenomena like the skin effect to acountercurrent in an arbitrary wire or winding [21]. These phenomena are already

3

1 Introduction

contained in the formulation of the full Maxwell equations.

In the case of transformer coils, a turn is often built with a thin, convoluted wirestructure, which leads to an almost homogeneous distribution of the current densityinside of the turn. The usual skin effect can be neglected, and induction currents donot occur. Coils without simulating the induction parts in the full Maxwell equationsare called stranded windings. In comparison, massive windings contain all inductionphenomena as in solid parts.

OutlineThe goal of this thesis is the integration of the introduced winding model inside ofthe full Maxwell system and furthermore the implementation and validation of theadjusted solver for the different winding types.

In the second chapter, the already existing mathematical formulation of the robust,full Maxwell system and its discretization is described and derived. Further, a shortoverview of the implemented solving techniques is depicted. Section 3 and 4 show therealization of the new winding model for both, stranded and massive windings. Themathematical inclusion inside of the full Maxwell system and also the implementationof the new structure is described exhaustively. The penultimate Section 5 deals withthe validation and comparison of the new model. It contains comparisons to theanalytical solution, refinement studies and multiple model comparisons of realisticbusiness simulations with other commercial solvers.

4

2 Robust Full Maxwell


2.1 Mathematical Formulation

The basic relations for all electromagnetic phenomena are the Maxwell equations. Fromthe general equations, a solvable numerical formulation is necessary, which is definedby [7] for the entire frequency domain. A short derivation of the used formulation isshown, which starts from the linear form of the Maxwell equation system.

∇ ·B = 0 (1a)

∇ · (εE) = ρ (1b)

∇× E = −∂B∂t

(1c)

∇× (1

µB)− ∂εE

∂t= j (1d)

With the linear material relations and variables

B = µ ·H and D = ε · E . (2)

• B[T ] is the magnetic flux

• H[Am

] is the magnetic field

• µ [Hm

] is the permeability

• ρ[ Cm3 ] is the charge density

• D[Asm2 ] is the electric flux

• E[ Vm

] is the electric field

• ε [ Fm

] is the permittivity

• j[ Am2 ] is the current density

For the solving procedure in the FEM-solver, the formulation above will be adjustedso that the resulting LES is well defined and conditioned for the entire frequencydomain. A scalar potential ϕi and a vector potential A are introduced for the solutionvariables [7]. ϕi symbolizes the induction component of the scalar potential, which canbe considered in this Section as ϕ.

B = ∇×A (3a)

E = −∇ϕ− ∂A

∂t(3b)

This equations are inserted in system (1) and further the Coulomb gauge condition isincluded. [7].

5


−∇ · (ε∇ϕ)− ∂

∂t∇ · (εA) = ρ (4a)

∇× (1

µ∇×A) + ε∇∂ϕ

∂t+ ε

∂2A

∂t2= j (4b)

∇ · (εA) = 0 (4c)

The four first-order partial differential equations from the known Maxwell formulationare transferred to two-second order equations, which simplifies the computation. Thementioned Coulomb equation is necessary to defined the potentials A and ϕ uniqueand solvable. [2]

For a given frequency an electromagnetic property has periodic oscillations so that forsolubility all parts of the equation system above are written in a frequency formulation.Because of the periodic behavior, every variable can be considered with a spatial anda transient part. [14] The oscillations also hold for the introduced potentials A and ϕso that the computation can be reduced to the spatial part . [2]

ρ = Re[ρ(x) · exp(iωt+ Φ)] (5a)

j = Re[j(x) · exp(iωt+ Φ)] (5b)

A = Re[A(x) · exp(iωt+ Φ)] (5c)

ϕ = Re[ϕ(x) · exp(iωt+ Φ)] (5d)

The Equation below shows the transformation from the solved spatial part back in thefrequency domain.

ϕ(X, t) =

∫ ∞−∞

ϕ(X, t) · exp(iωt) dω (6)

Combining the formulations (4) and (5) results in a system, where only the complexspatial parts A, ϕ, ρ and J remain.

∇× (µ−1 ∇× A) + iεω∇ϕ− εω2A = j (7a)

∇ · (εA) = 0 (7b)

−∇ · (ε∇ϕ) = ρ (7c)

6


In addition, the specific Ohm’s law (8) and the continuity equation (9) are consideredand assumed to be fulfilled [2].

j = σE = −σ∇ϕ− iσωA (8)

iωρ+∇ · j = 0 (9)

The index above all spatial values in the formulation is neglected and the variable ϕis due to a robust formulation in low frequency domain considered as ϕ+ ψ [7]. Withthe Equation (8) and Equation (9) the robust Maxwell formulation finally results to

∇× (1

µ∇×A)− (ω2ε− iωσ)A

+(iωε+ σ)∇(ϕ+ ψ) = 0 , (10a)

∇ · (εA) = 0 , (10b)

∇ · (ε∇(ϕ+ ψ)) = 0 . (10c)

This robust formulated system is written in weak form by Hiptmair [7].

〈µ−1∇×A,∇×A′〉 − 〈(ω2ε− iωσ)A,A′〉+〈(iωε+ σ)∇ϕ,A′〉+ iω〈ε∇ψ,A′〉 = 0 , (11a)

〈εA,∇ϕ′〉 = 0 , (11b)

〈ε∇ϕ,∇ψ′〉+ 〈ε∇ψ,∇ψ′〉 = 0 . (11c)

Solving the robust Maxwell system with a numerical method requires an explicit prob-lem formulation so that for the formulation (11) the following bilinear formulations areintroduced.

a1(A,A′) =1

µ

∫Ω

∇×A ∇×A′dx−∫

Ω

(ω2ε− iωσ)AA′dx (12a)

a2(ϕ,A′) =

∫Ω

(iωε+ σ)∇ϕA′dx (12b)

a3(ψ,A′) =

∫Ω

iωε∇ψA′dx (12c)

a4(A, ϕ′) =

∫Ω

εA∇ϕ′dx (12d)

a5(ϕ, ψ′) =

∫Ω

ε∇ϕ∇ψ′dx (12e)

7


a6(ψ, ψ′) =

∫Ω

ε∇ψ∇ψ′dx (12f)

Further, the test functions A, ϕ and ψ are defined in Sobolev spaces [7].

V := v ∈ H(curl,Ω) : ∇× vt = 0 on ∂Ω,

∫τ

v dS = 0 (13)

H(U) := ψ ∈ H1(Ω) : ψΓ0 = 0, ψΓ1 = U (14)

This finally leads to the variational formulation of the robust full Maxwell equationsystem.

Find A′ ∈ V , ϕ′ ∈ H(U) and ψ′ ∈ H1e (Ω) so that

a1(A,A′) + a2(ϕ,A′) + a3(ψ,A′) = 0 , (15a)

a4(A, ϕ′) = 0 , (15b)

a5(ϕ, ψ′) + a6(ψ, ψ′) = 0 . (15c)

2.2 Numerical Discretization

For the numerical solving of the variational formulation, a finite element method basedon the Galerkin method is used inside of the solver [6]. The geometry is discretized incurved tetrahedral of second order, where every element contains ten nodes - with fourmain nodes at the corners - and six edges. A visualization of a second order tetrahedralis depicted in Figure 3.

The two introduced potentials ϕ and A are defined on every tetrahedral. For the scalarpotential exist nodal, shape functions Wi on the main nodes, whereas first-order, edgefunctions Vi are used for the vector potential at the six edges [7]. A curved edgefunction can be described by the expression below [2].

VA B = WB · ∇WA −WA · ∇WB (16)

If the numbers of used edges m, interior nodes n and outer nodes o are known, thediscretization for the test functions can be defined with the coefficients αi, βi, γi andthe corresponding shape functions.

A =m∑i=1

αiVi , with Vi ∈ V (17a)

8


5

1

3

10

6

4

9

2

7

8

Figure 3: 3D Tetrahedral Element

ϕ =n∑j=1

βjWj , with Wi ∈ H(U) (17b)

ψ =o∑

k=1

γkWk , with Wi ∈ H1e (Ω) (17c)

Inserting the approaches (17) in the variational formulation (15) leads to the discretizedequation system of the full Maxwell equations.

m∑i=1

αi · a1(Vi,Vl) +n∑j=1

βj · a2(Wj,Vl)

+o∑

k=1

γk · a3(Wk,Vl) = 0 , l ∈ (1,m) (18a)

m∑i=1

αi · a4(Vi,Wl) = 0 , l ∈ (1, n) (18b)

n∑j=1

βj · a5(Wj,Wl) +o∑

k=1

γ · a6(Wk,Wl) = 0 , l ∈ (1, o) (18c)

This discretized equation system of formulation (11) is written in a solvable discrete

9


matrix formulation. A1 A2 A3

A4 0 00 A5 A6

·Aϕψ

=

000

(19)

2.3 Solving Procedure

The resulting matrix system Ax = b has to be solved with a numerical algorithm. Forthis several convenient solving methods are considered, whereas the choice depends onthe simulation settings and the model discretization.

Matrix AssemblingThe existing solver can simulate models with up to 14 million elements [12], which leadsto a huge and sparse system matrix, where only interacting nodes and edges cause amatrix entry. The memory consumption, while saving the matrix can be improved byusing the compressed sparse row data format, which consists of only three vectors. Thefirst vector holds all non-zero values of the matrix in a row-wise order but contains noinformation about the position. Hence the second vector saves the column index ofevery non-zero entry, whereas the last vector keeps the index of a row in the two othervectors. This new storage scheme provides a speed up in the matrix computations andneeds less storage. [15]

Solving with a Direct SolverA small and particularly regular system matrix can be solved with a direct solver,whereas inside of HADAPT the parallel direct sparse solver PARDISO is used [16].Indeed, in the case of solving system (18) the solving with a direct solver becomes dueto singularities in the first equation unstable so that the usage of an iterative solveris required [8]. Furthermore the direct solving needs more memory so that it is onlyusable for small or coarse meshed geometries.

Solving with the BiCGSTAB-MethodAs described in [8] an iterative solver has to be considered for singular matrices. Nu-merical iterative solving methods based on the Krylow-subspace save only the sparsematrix and a few computation vectors, which leads in general to a low memory con-sumption [5]. One of those Krylow methods, the conjugate gradient algorithm, leadsto an exact solution of the LES Ax = b after k iterations, where k is the dimensionnumber of the matrix A. The requirement for the usage of the primary CG-algorithm isa sparse and symmetric matrix, so that an approvement leads to the BiCG-algorithm,which also allows unsymmetric matrices, but exhibits bad stable conditions [17]. Thefurther improved stabilized biconjugate gradient algorithm (BiCGStab) by van der

10


Vorst stabilizes this Krylow method and is used for the more complex iterative solvingin HADAPT [19].

Solving with the Schur-ComplementDepending on the problem type and the applied boundary conditions the resultingequation system have a different form due to additional, extended equations for exci-tations. The extended dense matrix part increases the memory and deteriorates the sol-ubility so that a new solving technique is required. Besides, a separation of the solvingis in particular cases more efficient. Hence the sparse matrix part without excitationsis detached and solved in a previous solving step with the usual BiCGStab method.The extended rows and columns are solved with the particular Schur -algorithm, whichsolves the LES iteratively with a better efficiency [3]. The matrix and vectors of theequation system are divided into the following parts,

(M B1

B2 C

)·(xI

)=

(bU

), (20)

where the first row can be extracted.

M · x+B1 · I = b

⇔ x = M−1 · (b−B1 · I) (21)

If x is inserted into the second row, it follows to

B2 · x+ C · I = U

=⇒ B2 ·M−1 · (b−B1 · I) + C · I = U

=⇒ (C −B2 ·M−1 ·B1) · I = U −B2 ·M−1 · b (22)

Finally, the system is solved according to the procedure described in [8].

1. Compute with BiCGStab: br = M−1b

2. Compute with BiCGStab: bi = M−1B1i for every column i of the matrix B1.This leads to a new matrix Bi consisting of the solution vectors bi.

3. Compute L := (C−B2Bi)

4. Compute with direct PARDISO Interface: L−1

5. Solve with direct PARDISO Interface: I = L−1 (U−B2br)

6. Compute x = br −BiI

11


For every iterative solution method above a preconditioner is applied to reach a bettercondition number and convergence of the LES. Depending on the solving method thepreconditioner changes. [8] With a well-conditioned problem, the numerical algorithmof the solver converges in less than ten steps to residuals of 10−14.

12

3 Excitation of Stranded Windings


The basic equation system of the full Maxwell formulation (11) is considered, due to anintegration of the new winding model. In the case of a prescribed winding excitationthe current density in the winding changes regardless of the type of excitation. Becauseof the additional current density j0 in the winding Ohm’s law from the derivation inSection 2.1 have to be extended. It leads to

J = σE + j0 = −σ∇ϕi − iσωA + j0 . (23)

As described in Section 1 the distribution in a stranded winding is considered as ho-mogeneous. Therefore all induction effects are neglected, which can be expressed byan electric conductivity of zero (σ = 0). With this modification, the extended Ohm’slaw and some conversions the equation system for a stranded winding excitation endup with

∇× (µ−1 ∇×A)− ω2εA + iωε∇(ϕ+ ψ) = j0 , (24a)

∇ · (εA) = 0 , (24b)

∇ · (ε∇(ϕ+ ψ)) = 0 , (24c)

or in week formulation

〈µ−1∇×A,∇×A′〉 − 〈ω2εA,A′〉+〈iωε∇ϕ,A′〉+ iω〈ε∇ψ,A′〉 = 〈j0,A′〉 , (25a)

〈εA,∇ϕ′〉 = 0 , (25b)

〈ε∇ϕ,∇ψ′〉+ 〈ε∇ψ,∇ψ′〉 = 0 . (25c)

The further formulation of the equation system and its discretized form depend on theexcitation type of the winding. In the first step, the more intuitive current excitationwith a prescribed I0 is discussed. In the second section, the voltage excitation U andthe associated modifications are described.

The current density j0 leads to an analogous scalar potential ϕ0 in the winding domain.Hence, the entire potential ϕ end up with ϕ = ϕ0 + ϕi.

3.1 Current Excitation

In the case of a current excitation the prescribed value is the current I0, which consistsof a peak value I0 ∈ R and due to different phases an attendant phase angle Φ ∈ R.The conversions and effects to the equation system are described in the following.

13


3.1.1 Mathematical Formulation

Equation system (25) shows, that the influence of windings is an additional currentdensity in the right hand side. For a mathematical formulation of the current excitationthe relation between j(x) and the prescribed I0 is required. In general the current I(t)is defined by the complex and transient relation below.

I(t) = Re[I0 · eiωt] , I0 ∈ C (26)

I0 = I0 · eiΦ , I0,Φ ∈ R (27)

= I0, Re + i · I0, Im (28)

The current I(t) consists of the complex value I0, that depends on the phase angle Φand the prescribed peak value I0. Furthermore the transient term eiωt can be neglectedduring the computation as seen in Equation (5). For an arbitrary surface A the complexcurrent I0 ∈ C through this surface can be derived with a different definition.

I0 =

∫A

j0(x) dA , j(x, t) ∈ C (29a)

=

∫A

(j0, Re + i · j0, Im) dA (29b)

=

∫A

j0, Re dA+ i ·∫A

j0, Im dA (29c)

= j0, Re A+ i · j0, Im A (29d)

A comparison of both approaches results in a relation between the current density j0

and the prescribed current I0.

j0, Re =I0, Re

A=I0 · cos(Φ)

A(30a)

j0, Im =I0, Im

A=I0 · sin(Φ)

A(30b)

With equations (30) the complex current density j(x) can be calculated. As describedin the introduction a modeled cylinder can symbolize a coil with an arbitrary amountof turns. Hence the applied physical current of a usual coil has to be modified for asimplified cylinder. More precisely the applied current in the cylinder is the sum of allcurrents in the turns.

Icylinder = Nturns · Icoil (31)

14


With the weak formulation (25) the new variational formulation of the current excita-tion results in

Find A′ ∈ V , ϕ′ ∈ H(U) and ψ′ ∈ H1e (Ω) so that

a1(A,A′) + a2(ϕA′) + a3(ψ,A′) = f(A′) , (32a)

a4(A, ϕ′) = 0 , (32b)

a5(ϕ, ψ′) + a6(ψ, ψ′) = 0 . (32c)

The additional bilinear form of the right hand side f(A′) contains the current density,which can be expressed by the current as shown in Equation (30). The direction ofthe current density is due to a prescribed excitation always tangential in direction eϕ.Hence, a vector current density j can be derivated.

f(A′) =

∫Ω

j0(x) ·A′ dx , (33)

=

∫Ω

I0

A· eϕ(x) ·A′ dx . (34)

In the following the discretization and quadrature of the new Equation (34) are depictedmore precisely.

3.1.2 Numerical Discretization

Based on the variational formulation (32) the new discretized formulation yields to

m∑i=1


βj · a2(Wj,Vl)

+o∑

k=1

γk · a3(Wk,Vl) = f(Vl) , l ∈ (1,m) , (35a)

m∑i=1

αi · a4(Vi,Wl) = 0 , l ∈ (1, n) , (35b)

n∑j=1


k=1

γ · a6(Wk,Wl) = 0 , l ∈ (1, o) . (35c)

This leads to a new matrix formulation, where the adjustment is the additional right

15


hand side term A1 A2 A3

A4 0 00 A5 A6

·Aϕψ

=

f00

(36)

with

f =

f(V1)...

f(Vm)

∈ C(m×1) .

For every edge in the winding domain, the bilinear form f(Vi) have to be computed.Hence, a four-point quadrature method is used on every tetrahedral to calculate theentry of the right hand side.

f(Vi) =

∫Ω

j0(x) ·Vi dx ∀i ∈ (1, NEdges) (37a)

≈NTets∑j=1

∫Tetj

j0(x) ·Vi dx ∀i ∈ (1, NEdges) (37b)

The integral over every tetrahedral is computed with a quadrature method as describedin [13]. ∫

Tetj

j0(x) ·Vi dx =

∑Nqp

k=1 ωk · j0(xk) ·Vi Edgei ∈ Tetj0 Edgei 6∈ Tetj

(38)

This leads to the final equation for the quadrature.

f(Vi) ≈NTets∑j=1

Edgei∈Tetj

Nqp∑k=1

ωk · j0(xk) ·Vi ∀i ∈ (1, NEdges) (39)

With a calculated j0(x) from the prescribed I0, the quadrature can be computed andfilled in the right hand side vector to solve the LES.

16


3.1.3 Implementation

Equation (39) shows that the quadrature for every edge i in the winding domain is asum over all tetrahedral, which contain this edge i. This would lead to a complexity ofO(No_Edges× No_Tets) as depicted in the pseudo code in Figure 4a. Instead of usingthis procedure the complexity can be reduced by iterating with an outer loop over alltetrahedral and calculate all edge quadratures of an element at the same time. Thisimproves the complexity to O(No_Tets), which is depicted in the code in Figure 4b.After the calculation, the respective value will be added in the right hand side.

1 f o r i = 1 : No Edges2 i f edgeIsNotInWinding cont inue ;3 rhs pos = getIndex ( i ) ;45 f o r j = 1 : No Tets6 i f edgeIsNotInTet cont inue ;7 rhs [ rhs pos ] += getBLF ( Tet j ) ;

(a) code with outer edge loop

1 f o r i = 1 : No Tets2 i f tetIsNotInWinding cont inue ;3 j t e t = getAllBLF ( Tet i ) ;45 f o r j = 1 :66 rhs pos = getIndex ( Tet i , j ) ;7 rhs [ rhs pos ] += j t e t [ j ] ;

(b) code with outer tet loop

Figure 4: Code for Right Hand Side Calculation

3.2 Voltage Excitation

The second excitation type is a prescribed voltage drop U ∈ R with the equivalentphase angle Φ ∈ R. The similar required relation between the voltage U and thecurrent density j0 is, in this case, more elaborate and discussed in the following Section.

3.2.1 Mathematical Formulation

The relation between the prescribed voltage and the required current density j0 canbe derived by separating the total prescribed voltage U . Besides the applied windingvoltage U0, there exists an additional induction component Ui.

U = U0 + Ui =L

A · σ· I0 +

L

Vring

∫Ω

E · eϕ dV (40a)

=L

A · σ· I0 +

L

Vring

∫Ω

(−∇ϕi − iωA) · eϕ dV (40b)

=L

A · σ· I0 −

L

Vring

∫Ω

iωA · eϕ dV (40c)

The winding voltage U0 can be defined with the knowledge of the resistance R, whereL is the perimeter, A the cross-section and σ the electric conductivity.

17


R =L

A · σ

The second term in Equation (40) is the induced voltage Ui. As seen the inducedvoltage is an integration over the electric field, which consists of two components. Thefirst component ∇ϕi is due to the circular structure of a winding always zero andtherefore neglected.

Proof ∫Ω

∇ϕi · eϕ dV (41a)

=

∫ Γ2

Γ1

∇ϕi · eϕ dS Γ1 = Γ2 = A (41b)

= ϕi(Γ2) · eϕ(Γ2)− ϕi(Γ1) · eϕ(Γ1) (41c)

= 0 2 (41d)

The remaining second part Ui leads together with the applied voltage U0 to the totalvoltage U .

U = U0 + Ui = R · I0 −L

Vring

∫Ω

iωA · eϕ dV (42)

The voltage drop of a spiral coil Ucoil has to be adjusted analogous to the current inEquation (31). Because of a shortened length of a simplified cylinder, the prescribedvoltage has to be divided by the number of turns.

Ucylinder =Ucoil

Nturns

(43)

The relation between U and I0 leads to an additional equation for every winding in thefull Maxwell equation system. Therefore Equation (42) is expressed by new introducedbilinear formulations and the relation B2 · x + C · I0 = U . Further, the expression ofthe current density B1 in the right hand side is transformed with an additional bilinearform to the matrix side because the required current I0 is now a solving variable in thesolution vector and no longer prescribed.

18


B1(I0,A′) = −

∫Ω

eϕA·A′ dx (44a)

B2(A, I ′0) = − L

Vring

∫Ω

iωA · eϕ dx (44b)

C(I0, I′0) =

L

AR · σ(44c)

With the formulation (15) and the additional bilinear formulations the extended vari-ational formulation of the equation system with voltage excitation results in

Find A′ ∈ V , ϕ′ ∈ H(U) andψ′ ∈ H1e (Ω) so that

a1(A,A′) + a2(ϕ,A′) + a3(ψ,A′)−B1(I0,A′) = 0 , (45a)

a4(A, ϕ′) = 0 , (45b)

a5(ϕ, ψ′) + a6(ψ, ψ′) = 0 , (45c)

B2(A, I ′0) + C(I0, I′0) = U . (45d)

3.2.2 Numerical Discretization

The variational formulation (45) is again transfered to a discretized equation system,whereas Vi and Wi are the shape functions from Equation (17). Further Nvoltage

expresses the number of windings with stranded voltage excitation.

m∑i=1


βj · a2(Wj,Vl)

+o∑

k=1

γk · a3(Wk,Vl) +

Nvoltage∑h=1

I0hB1(I0h ,Vl) = 0 , l ∈ (1,m) , (46a)

m∑i=1

αi · a4(Vi,Wl) = 0 , l ∈ (1, n) , (46b)

n∑j=1


k=1

γ · a6(Wk,Wl) = 0 , l ∈ (1, o) , (46c)

m∑i=1

αi ·B2(Vi, I0l) + C(I0l , I0l) = Ul , l ∈ (1, Nvoltage) . (46d)

This discretized equations lead to a global matrix system, which describes (45) in a

19


discrete solvable formulation.A1 A2 A3 B1

A4 0 0 00 A5 A6 0B2 0 0 C

·AϕψI0

=

000U

(47)

3.2.3 Implementation

The complex-valued entries of the additional matrices B1, B2 and C are calculatedagain with the usage of four-point quadrature methods. Equation (45) shows, that theimaginary part of B1 and C and the real part of B2 are zero and therefore neglected inthe quadratures. The small matrix C results in a diagonal matrix, whereby the entriesof matrix C are computed directly without any quadrature.

Both, Matrix B1 and B2 consists of entries for which the calculation of the bilinear form∫Ωeϕ ·A dx is required. This is an analogues quadrature to the method of

∫Ωj0 ·A dx

in Section 3.1 and can be computed with the knowledge of the edge function and thewinding tangential eϕ at every quadrature point. Because of complexity reasons thequadrature is performed in an outer loop for all tetrahedral. The pseudo code of thefour-point quadrature is visualized in Figure 5.

1 f o r edge j = 1 :62 e dg e fu nc t i o n s = getEdgeFunctionsAtQps ( t e t i , edge j ) ;34 f o r k = 1 : No Quadrature Points5 e ph i = getEPhi ( p o s i t i o n ( k ) ) ;6 edge = ed g e fu nc t i o n s [ k ] ;78 r e s u l t [ j ] += weight [ k ] ∗ det [ k ] ∗ edge ∗ e ph i ;

Figure 5: Code for Quadrature∫Teti

eϕ ·A dx

The calculation of the winding tangential eϕ is described in detail in the next Section.

3.3 Calculation of Winding Tangential eϕ

As seen either for the functional f(V) or the bilinear forms B1(I0,V) and B2(V, I0)the knowledge of the winding tangential is necessary. In a modeled stranded windingthe current is enforced in a circular, tangential direction, which is homogeneous alongan arbitrary cross-section. The following two sections describe the analytical approachand a numerical solution for the winding tangential eϕ.

20


3.3.1 Failure of Analytical eϕ

The analytical eϕ points always in the tangential direction of the winding. With theknowledge of geometric parameters, the winding tangential can be calculated at anarbitrary point. With the discrete mesh data, the midpoint of the winding ~xwinding,the midpoint of an arbitrary cross-section ~xarea and the normal vector ~n to this areaare simply calculable. Further, the winding tangential eϕ is computed at an arbitrarywinding point ~p.

~rref = ~xwinding − ~xarea (48a)

~a = ~rref × ~n (48b)

~r = ~p− ~xwinding (48c)

~rproj = ~rref ·~r · ~rref

||~rref||2+ ~n · ~r · ~n

||~n||2(48d)

eϕ = ~rproj × ~a (48e)

The calculated vector ~a describes the orthogonal winding axis, cutting the midpointalong the symmetry line. Hence, the resulting winding tangential eϕ depends on thevector ~a and the projected radius ~rproj, and always lies in the winding plane.

Although the winding tangential have a geometrical approach, the equation systemleads to a not solvable right hand side. In system (24) the current density and thereforethe winding tangential eϕ influence the new derived right hand side in the first equation.If the divergence is multiplied to this equation, it results in

∇ · (∇× (µ−1 ∇×A))−∇ · (ω2εA) +∇ · (iωε∇(ϕ+ ψ)) = ∇ · j0 . (49)

For the observed case of stranded windings, the conductivity σ is modeled as zero.With ω = 0 and a relative permeability of µ = 1 the left side of the equation aboveremains with ∇ · (∇× (∇×A)), which end up with zero. For the proof the potentialA = ∇×A is introduced.

Proof

∇ · (∇× (∇×A)) = ∇ · (∇× A) = ∇×

∂Az

∂y− ∂Ay

∂z∂Ax

∂z− ∂Az

∂x∂Ay

∂x− ∂Ax

∂y

(50a)

=∂2Az∂x∂y

− ∂2Ay∂x∂z

+∂2Ax∂y∂z

− ∂2Az∂x∂y

+∂2Ay∂x∂z

− ∂2Ax∂y∂z

(50b)

21


= 0 2 (50c)

With this proof the left side of Equation (49) is zero so that also ∇·j0 = j0 ·∇·(eϕ) hasto be zero. If the winding tangential eϕ is expressed in cylindric coordinates eϕ(r,Φ, z)the divergence evaluation can be simplified [10].

∇ · (eϕ) =1

r(eϕr)r +

1

reϕΦ,Φ

+ eϕz,z = 0 (51)

All derivatives of the components above are zero so that the analytical winding tan-gential eϕ and further the first equation of system (24) are divergence-free. Instead ofsolving the analytical equation a discretization is applied for using the finite elementmethod. With this discretization, the divergence of the derived winding tangential inEquation (48) is not zero, which leads to a not solvable right hand side and problemsin the solving procedure. For a better solvability of the equation system, the wind-ing tangential is computed numerically, which results in a numerical, divergence-freewinding tangential eϕ.

3.3.2 Discrete divergence-free eϕ

As shown the winding tangential has to be determined numerically, so that ∇ · eϕ = 0holds also for the numerical discretization. A new ansatz introduces a scalar potentialϕ0 and an additional relation between the winding tangential eϕ and the potential ϕ0.For this derivation, the gradient of the unknown potential ϕ0 is observed in cylindriccoordinates.

∇ϕ0 =∂ϕ0

∂rer +

1

r

∂ϕ0

∂ΘeΘ +

∂ϕ0

∂zez (52)

Due to a wanted constant tangential eϕ Equation (52) should end up with eΘ. Thisleads to a potential ϕ0 = Θ and an additional scaling factor r so that the approach forthe winding tangential eϕ end up with

eϕ = r · ∇ϕ0 . (53)

The used radius r and the potential ϕ0 = Θ ∈ [0, 2π] lead to a constant distributionof the winding tangential. Appling this approach causes the same numerical problemsbecause eϕ have to be divergence-free on a discrete level. Therefore the followingequation has to be solved which ensures a divergence-free eϕ.

⇒ ∇ · (eϕ0) = ∇ · (r · ∇ϕ0) = r ·∆ϕ0!

= 0 . (54)

22


In the following section, the introduced winding potential ϕ0 is written as ϕ, becauseother potentials can be neglected in the solving procedure of ϕ0.

Mathematical FormulationThe formulation (54) is a well-posed Poisson problem, which can be solved with Neu-mann boundary conditions. Because of the geometry of a winding, a discontinuityof the potential ϕ exists at a defined cross-section A, where the Neumann boundaryconditions of the problem can be defined. The area, depicted in Figure 6, is divided intwo surfaces, Γ1 and Γ2, with different applied boundary conditions.

Γ1

Γ2

ϕ2ϕ1

Figure 6: Abstract Cross-Section of a Winding

For a unique solution of the potential ϕ one point in the entire system is defined as aDirichlet point with ϕ = 0, which is considered on the surface ΓD. (ΓD ⊂ Γ2)

The local, linear testfunction ϕ′ is defined at all mainnodes by

ϕ′ ∈ H1(Ω) with ϕ′ = 0 on ΓD ,

which leads to the equation ∫Ω

r ·∆ϕ · ϕ′ dV = 0 (55)

or with partial integration∫

Ωϕ · ∇ · (~v) dV =

∫∂Ω~n · ϕ · ~v dS −

∫Ω~v · ∇ϕ dV to

∫Ω

∇ · (r · ∇ϕ) · ϕ′ dV =

∫δΩ

~n · r · ∇ϕ · ϕ′ dS −∫

Ω

r · ∇ϕ · ∇ϕ′ dV . (56)

23


With Equation (55) the formulation of the Poisson problem end up with∫Ω

r · ∇ϕ · ∇ϕ′ dV =

∫δΩ

~n · r · ∇ϕ · ϕ′ dS , (57)

~n · r · ∇ϕ = ~n · eϕ =

1 on Γ1

−1 on Γ2\ΓD0 on Ω\(Γ1 ∪ Γ2)

(58)

The dirichlet boundary condition on ΓD can be negleted, because ϕ′ = 0 on ΓD. Withthis relations the introduced bilinear and variational formulations are

a(ϕ, ϕ′) =

∫Ω

r · ∇ϕ · ∇ϕ′ dV , (59a)

f(ϕ′) =

∫Γ1

ϕ′ dS −∫

Γ2

ϕ′ dS . (59b)

⇒ Find ϕ′ ∈ H1(Ω) so that

a(ϕ, ϕ′) = f(ϕ′) . (59c)

Numerical DiscretizationWith the definition of first-order shape functions Φi and the number of winding nodesN the following discrete approaches result.

ϕ =N∑i=1

αi · Φi , (60a)

ϕ′ = Φi , i ∈ (1, N) . (60b)

Besides, the equation system is formed by

m∑i=1

αi · a(Φi,Φl) = f(Φl) , l ∈ (1,m) , (61)

with the attendant LESa(Φ1,Φ1) · · · a(Φm,Φ1)...

. . ....

a(Φ1,Φm) · · · a(Φm,Φm)

·α1

...αm

=

f(Φ1)...

f(Φm)

(62)

24


ImplementationDuring the implementation, the detection and indexing of the cross-section with therespective boundary conditions is the most difficult task. A closer look at the dis-cretization of the cross-section is depicted in Figure 7.

1

2

3

4

5

7

8

9

10

6

1

2

3

4

5

11

12

13

Γ1 (domain 1) Γ2 (domain 2)

Figure 7: 2D Discretized Cross-SectionNote: Domain Nodes, Neumann Boundary Nodes, Dirichlet Boundary Node

A physical mesh node on the surface belongs to both abstract surfaces and becomestwo nodes in the LES with different, opposite boundary conditions. All Neumannconditions are depicted yellow, and the one Dirichlet node with ϕ = 0 is visualizedorange. For better recognition, the nodes on the surface, as well as the winding nodes,are marked in arrays at the beginning. Further, the node interactions have to bestored in a two-dimensional array because only adjacent nodes are taken into accountin the further calculation, whereas there is no interaction between the discontinuity inFigure 7.

It is essential that the additional nodes on Γ2 obtain due to a better distinction anindex at the end of the usual node array. The assignment to the surfaces can beproceeded in practice by distinguishing the material name.

The computation of all matrix entries requires the knowledge of the projected radius~rproj at every quadrature point. This radius is already calculated in Equation (48) andcan be used in the quadrature method.

~rproj = ~rref ·~r · ~rref||rref ||2

+ ~nref ·~r · ~nref||nref ||2

(63)

With this radius, the quadrature method is able to compute the bilinear form a(Φi,Φl),before filling the CRS matrix with the results. The resulting matrix M is symmetric,because of the symmetry of a(Φi,Φj) and has due to the few node interactions a sparse

25


form.

For the right hand side vector, a surface integral over (Γ1∪Γ2)\ΓD is necessary. Insteadof an integral over entire elements, the calculated integral is only a surface integral,which requires a few adjustments.

If both sides of the symmetric equation system are defined, the algorithm PARDISOsolves the matrix-vector system due to a symmetric, small matrix within a directmethod as described in Section 2.3. The solution array ~α contains the solution for ϕat every node in the winding domain so that with the local gradients and the radius~rproj the winding tangential eϕ can be revealed by

eϕTeti=

∑MainNodej∈Teti

r · ϕ · localGradHatj . (64)

The calculated winding tangential eϕ results in a numerical divergence-free variation.Inserting the pre-calculated eϕ in the quadratures of the global system from Section 3.1and 3.2 leads to a solvable full Maxwell equation system, which is solved subsequently.The full process is visualized in the flowchart in Figure 8.

3.4 Postprocessing Calculations

After solving the equation system of the full Maxwell formulation and all field val-ues, postprocessing values like the induced voltage, ohmic losses or acting forces arecalculated.

Induced VoltageIn the case of a prescribed voltage excitation, the voltage drop in a winding is al-ready known, while for prescribing a current the induced voltage is computed withEquation (40) in a postprocessing step.

U0 =L

AR · σ· I0 (65a)

Ui = − L

Vring

∫Ω

iωA · eϕ dV (65b)

For the applied voltage all values are prescribed. The second, induced componentcontains the solving potential A, which is in a post-processing step available. The sumof both voltage components calculates the total winding voltage U = U0 + Ui.

26


Read mesh data, material parameter andboundary conditions from the input file

Detecting winding struture, assemblingof matrices and calculation of entries

Solving equation system withthe direct PARDISO method

Saving solution in global array phi forfurther winding tangential calculation

Assembling of full Maxwell system, calcu-lation of matrix and right hand side entries

Solving equation system with the Schurcomplement technique and BICGStab

Saving all solved potentials in global arrays

Calculation and visualization of all physical field valuesand postprocessing properties (losses, forces, voltages)

Figure 8: Flow Chart of Implemented Structure in HADAPTNote: HADAPT, Winding System, FullMaxwell System

Ohmic LossesFurther the losses in a winding are required for electromagnetic simulations, whereby

ohmic losses are defined by the expression P =j20·Vring

σ. The occurring current density

can be computed from the current, which is either prescribed or solved in the equationsystem. The electric conductivity is also given, and the domain volume is calculablewith the discretized geometry data. Hence the ohmic losses are computable in the entirewinding domain and can be used further as heat sources in additional or coupled CFD-simulations. In general ohmic losses occur in the form of heat and therefore influencethe fluid around. A one-way CFD coupling is processed with the commercial softwareFLUENT and solves the heat equation for the temperature distribution, based on theohmic losses. The temperature dependency of the electric conductivity could also be

27


considered in a two-way coupling process.

ForcesIn the last post-processing step the acting forces, which are induced by the currentdensity j0 are evaluated. The forces are based on the formulation j0×B and calculatedwith a quadrature method for every element. The magnetic field B is determined withthe relation ∇×A from definition (3), whereas the current density j0 is again eitherprescribed or solved in the full Maxwell equation system. The sum of all elementalcross products results in the acting force F.

28

4 Excitation of Massive Windings


The previous Chapter described the simulation of stranded windings, where all in-duction phenomena and therefore the equation parts with electric conductivity areneglected. Instead of this, massive windings depict all electromagnetic effects by mod-eling the entire equation system (11). Therefore the additional induction parts areconsidered in the following.

As seen in Section 3, the general equation system is extended by a winding currentdensity, which can be derived from Ohm’s law (23). Hence, the equation system formassive windings has an excitation term in the right hand side and additional inductionterms in the left matrix.

〈µ−1∇×A,∇×A′〉 − 〈(ω2ε− iωσ)A,A′〉+〈(iωε+ σ)∇ϕ,A′〉+ iω〈ε∇ψ,A′〉 = 〈j0,A′〉 , (66a)

〈εA,∇ϕ′〉 = 0 , (66b)

〈ε∇ϕ,∇ψ′〉+ 〈ε∇ψ,∇ψ′〉 = 0 . (66c)

The solving of a previous equation system for stranded windings is necessary becausethe used winding tangential has to be numerical divergence-free as described in Sec-tion 3.3.1. For massive windings, the first equation of system (66) is extended withthe induction terms which compensate numerical deviations of the winding tangential.Hence the formulation is stabilized, and the equation system should converge in thecase of applying an analytical winding tangential. This allows the usage of the analyt-ical winding potential ϕ0 without solving the previous equation system. Instead, ϕ0 iscalculated with the angle Θ in the winding plane.

As opposed to calculating a unit winding tangential as in Section 3 the expressionj0 = σ · ∇ϕ0 is used in a physical context. For massive windings, the value of σ is theactual winding conductivity. Besides, the resulting potential has a discontinuity withthe value of the winding voltage U . This simplifies the relation between voltage andcurrent and leads with a constant σ to an inhomogeneous current density distributionin the winding. Because of the induction effects, this distribution also holds for thecurrent density in a physical context. With the knowledge of the current density, theequation system for massive windings can be solved analog to the methods describedin Section 2.3.

4.1 Current Excitation

Due to the induction effects, the total current for current excitation is the sum of anapplied winding current I0 and the induced term Ii. This total current I is defined by

29


the following integral over an arbitrary cross-section of the winding.

I =

∫A

j dA (67a)

=

∫A

j0 + ji dA (67b)

=

∫A

−σ∇ϕ0 dA+

∫A

(−σ∇ϕi − σiωA) dA (67c)

= I0 + Ii (67d)

The applied current density j0 results from the defined potential ϕ0, whereby thepotential is considered in the following as a unit potential in the range of [0, 1]. Togetherwith the electric conductivity, the applied current density leads to

j0 = −σ · ∇ϕ0 = −σ · 1

2πr. (68)

ϕ2 = 1ϕ1 = 0 [ϕ] = 1

(a) unscaled unit potential ϕ0

ϕ2 = cϕ1 = 0 [ϕ] = c

(b) scaled unit potential c · ϕ0

Figure 9: Potential ϕ for Massive Voltage Excitation

Due to the unit potential ϕ0, which is depicted in Figure 9a, a scaling factor c ∈ C isintroduced and applied to all variables so that the resulting current I end up with theprescribed value of Iwanted. This scaling factor is defined by

c =Iwanted

I(69)

The factor c is used for scaling all occuring potentials.

ϕ0, scaled = c · ϕ0 (70a)

ϕi, scaled = c · ϕi (70b)

Ascaled = c ·A (70c)

With the scaled solving potentials the total current Iscaled end up with the prescribed

30


current Iwanted.

Iscaled =

∫A

−σ · ∇ϕ0, scaled dA+

∫A

(−σ · ∇ϕi, scaled − σ · iω ·Ascaled) dA (71a)

= c ·∫A

−σ · ∇ϕ0 dA+ c ·∫A

(−σ · ∇ϕi − σ · iω ·A) dA (71b)

With two introduced definitions

I0 =

∫A

−σ · ∇ϕ0 dA and (72a)

B · x =

∫A

(−σ · ∇ϕi − σ · iω ·A) dA (72b)

a new equation can be derived.

c · I0 + c ·B · x = Iscaled (73)

This current Iscaled is equal to the prescribed current Iwanted by definiton. Further, theformer matrix system Mx = b from Equation (66) is multiplied with the scaling factorc ∈ C.

c ·M · x = c · b (74a)

M ·(c ·Ac · ϕi

)= c · 〈j0,A〉 (74b)

Combining both equations, (71) and (74), yields to an extended system, which isassembled and solved inside of the solver.(

M −bB I0

)(c · xc

)=

(0

Iwanted

)(75)

Solving this equation system leads to the scaled variables of c · x = Ascaled, ϕi, scaledand the scaling factor c, which represents due to a unit potential ϕ0 ∈ [0, 1] thediscontinuity of ϕ0, scaled and further the voltage drop U . Figure 9b depicts the scaledpotential with a jump of the scaling factor c and therefore the voltage U .

Hence the induced voltage of a winding can be calculated with the discontinuity ofthe potential ϕscaled, whereby the part of ϕi, scaled can be neglected due to a constant

31


distribution along the winding perimeter.

U = c · [ϕ] = c · [ϕ0] = c (76)

With the extended equation system, the prescribed current value is set on the righthand side. The original equation system is transformed and scaled by embeddingan additional equation. The procedure of the discretization and implementation iscomparable to the process of stranded windings with voltage excitation in Section 3.2.

4.2 Voltage Excitation

If the prescribed value is the voltage U the mathematical formulation is simplified, dueto a neglected induction component Ui. Therefore the remaining U0 = U representsthe discontinuity of the potential ϕ0 and further the scaling factor c of the equationsystem. With Equation (74) the scaling factor c = U = U0 can be inserted in the righthand side without adding additional equations. The resulting solution potentials arealready scaled in a physical context.

M ·(Ascaled

ϕi, scaled

)= c · 〈j0,A〉 = U0 · 〈−σ · ∇ϕ0,A〉 (77)

This direct relation between the current density j0 and the prescribed voltage drop U ,leads to a simple calculation of the right hand side of matrix system (66). The proce-dure in the implementation is analog to the current excitation for stranded windings,where only a right hand side entry influences the full Maxwell system, without anyadditional equations.

4.3 Postprocesssing Calculations

Analog to Section 3.4 the values for the induced voltage, the ohmic losses and actingforces have to be calculated in a postprocesssing step.

Induced VoltageAs seen in Section 4.1 the winding voltage U can be expressed for massive windings bythe discontinuity of the potential ϕ0, sclaed. In the case of a current excitation, the jumpis represented by the complex factor c, which is solved inside of the equation systemand represents the voltage of the winding. So the total voltage U can be determinedby solving the extended equation system.

Ohmic LossesThe ohmic losses are computed analog to the case of stranded windings with the

32


relation P =j2·Vring

σ, whereby the difference for massive windings is the inhomogeneous

distribution of the current density j. Therefore a quadrature method have to be appliedwhere beside the winding current density j0 also the induced term ji = −σ(∇ϕi+ iωA)is considered. With this knowledge the omic losses in every element and further in theentire winding domain can be calculated.

ForcesThe last postprocesssing step calculates the forces, induced by the term j × B. Thedifference for massive windings is the additional induction current density ji, whichhave to be considered in the quadrature method. Hence the total current density endup with j = j0 − σ(∇ϕi + iωA), which leads with the magnetic field B to the actingforce.

Massive windings can be used in all models with different excitations and frequencies.Also, the combination of stranded voltage and massive current excitation inside ofone model is possible, which leads to a twofold extended equation system. Due to aninhomogeneous current distribution inside of a massive winding the simplification ofseveral turns to one cylinder cannot be performed in a physical context. Hence thenumber of turns is enforced to one for all massive windings.

In Section 5 the differences between massive and stranded windings are observed inmore detail in the entire frequency domain.

33

5 Validation

5 Validation

The implementation of the new winding model requires a validation, which is dividedand performed in three different steps. In a first step, the previous equation system forthe winding tangential of stranded windings is verified with a comparison to the ana-lytical solution. Subsequently, the entire winding model is validated for all excitationtypes and different geometries. Therefore some basic internal comparisons are per-formed with the already existing solver version. In the last step, a complex geometrywith practical benefit is modeled and validated with the commercial software MagNet.

5.1 Verification of the Winding Tangential Calculation

For the validation of the previous, Poisson equation system, solving the winding tan-gential eϕ for stranded windings, the analytical solution of the current direction hasto be derived. This solution can be used afterward for a comparison to the numericalresults.

5.1.1 Analytical Solution

The analytical approach for the previous equation system from Section 3.3.2 is

eϕ = r · ∇ϕ0 , (78)

whereas r is the projected radius in the winding plane and ϕ0 the introduced windingpotential. For a constant value of the variable r the resulting winding tangentialwould be higher on the interior wall of the winding, which is not desirable in the caseof stranded windings. In order to get a constant distribution the value for r is setto the projected radius in the winding plane. Further the relations ϕ0 = C · Θ andsubsequently ∇ϕ0 = ∆Θ

∆x= C·Θ

r·Θ eϕ are used for the derivation, with eϕ a unit vector inthe tangential direction.

eϕ = r · ∇ϕ0 (79)

= r · Creϕ (80)

= C · eϕ (81)

⇒ C = 1 (82)

So the analytical solution of the potential ϕ0 = C · Θ = Θ ∈ [0, 2π] have a linearprogression from zero until two pi. With a constant of C = 1 the analytical solution

34

5 Validation

of the winding tangential leads to the wanted, always tangential current direction eϕ.

eϕ = 1 · eϕ ⇒ |eϕ| = 1 (83a)

5.1.2 Compare Analytical and Numerical Solution of ϕ0 and eϕ

The analytical solutions of the potential ϕ0 and the winding tangential eϕ are derivatedin Section 5.1.1 and further used for a comparison to the numerical solution.

Comparison between Analytical and Numerical ϕ0

The analytical potential ϕ0 = Θ can be compared to the numerically computed solutionϕN0 from the previous, Poisson equation system for stranded windings. The results ofa winding with quadratic cross-section are depicted in Figure 10. The inner radius ofthe winding is 5 cm with a thickness of 1 cm.

(a) numerical solution ϕN0

N∑N

l=1 ||ϕ0l−ϕN

0l||∑N

l=1 ||ϕ0l||

353 0.010011718 0.002173758 0.0010016156 0.0009732757 0.00045128320 0.00033260637 0.00015

(b) relative euclidean error δϕ0

Note: The error is calculated in the L2-norm

Figure 10: Results of the winding potential ϕ0

Graphic 10a shows that the distribution of the computed ϕN0 is circular from zero totwo pi and therefore almost equal to the angle Θ. For a better comparison to theanalytical solution the relative error δϕ0 is calculated in Table 10b, where N is thenumber of winding tetrahedral.

The depicted relative error for different meshes shows that a mesh refinement resultsin a convergence of the numerical solution. Hence the small deviations in the relativeerror occur due to a discrete mesh structure.

Comparison between Analytical and Numerical eϕThe aim of solving the previous Poisson problem is the calculation of a numerical

35

5 Validation

divergence-free winding tangential eϕ. Based on the solution ϕ0 the numerical windingdirection can be computed by eϕ = r · ∇ϕ0 and is analyzed in Figure 11.

(a) numerical solution eNϕ

N∑N

l=1 ||eϕl−eNϕl

||2·Vl∑Nl=1 ||eϕl

||2·Vl

353 6.7491718 3.4713758 2.65216156 1.98532757 1.424128320 0.971260637 0.759

(b) relative euclidean error δeϕNote: The error is calculated in the L2-norm

Figure 11: Results of the winding tangential eϕ

The visualization of the numerical winding tangential eNϕ shows that the numericalsolution has variances around the desired magnitude of one. The differences can beobserved in more detail by calculating the relative error in Table 10b. Due to a vectorialeϕ the relative error δeϕ is calculated in the L2-norm.

For a refinement of the mesh, the relative error converges again to the analyticalsolution. This convergence becomes more visible in the convergence plot in Figure 12.

10 15 20 25 30 35 40 45 50 55 60

0.01

0.03

0.05

0.07c1 x-1

e

c2 x-1

Figure 12: Relative Error in L2 Norm in Log-log Plot of eϕNote: The two parallel lines with gradient of one determine the order of the convergence

The plots depict the L2-Norm of the error in a log-log scale for the shown winding

36

5 Validation

with different meshes. Therefore the relative error δeϕ is plotted against 3√N so that

a conclusion to the order of convergence is possible. The data points result in a linewith the gradient of minus one, which leads to an order of one for the convergence ofthe winding tangential eϕ.

Hence the previous equation system leads to good accordance and expectable results.

5.2 Internal Comparison With Existing Solver

After the validation of the previous equation system, the results of the entire fullMaxwell solver HADAPT are examined. In a first step, a comparison with the solveritself is accomplished. Therefore some pure windings with different setups are observed,and in addition, a comparison of two similar winding models with an additional wallis performed.

5.2.1 Winding Model With Different Simulation Settings

In this validation, a model with four independent winding turns is observed with dif-ferent simulation setups. The windings have arbitrary positions, peak values, phaseangles and are modeled either as stranded or massive turn, whereby Table 1 shows thedifferent setups precisely.

Excitation N I0 Φ

1 Stranded Current 1 100 A 0

2 Stranded Voltage 1 0.010427 V 74.65

3 Massive Current 1 100 A 0

4 Massive Voltage 1 0.009986 V 72.85

Table 1: Simulation Setup for Internal ComparisonNote: N is the number of turns in the winding

All windings have the same radius and a quadratic cross-section but are located atdifferent positions. The applied material is in all four cases copper with a relativepermeability and permittivity of one and an electric conductivity of 5.96 · 107 S

m. The

excitation of the four windings covers all different excitation types because two wind-ings are modeled as stranded, whereas the remaining are considered with the massivemodel. Besides, both voltage and current excitation are applied. The prescribed ex-citation of the voltage is computed in a previous simulation so that they correspondto an analog applied current of 100 A. Therefore the range of the current density inall windings should be comparable, whereas the distribution in the massive windingsvaries and presents a progression along the cross-section. With the relation j = I

Afor

stranded windings and an area A = 4 cm2 the current density can be computed. The

37

5 Validation

simulation is executed with a frequency of 50 Hz and a discretization of 800000 tetra-hedral elements. The following Figure shows the current density of all four windings.

(3)

(4)

(1)

(2)

Figure 13: Current Density j of Internal ComparisonNote: The number symbolizes the winding from Table 1

Figure 14a shows in more detail the results of a stranded winding. The depicted RMSvalue of the winding current density is in the range of 1.75 · 105 A

m2 , which can bederived analytical by calculating the RMS (Root Mean Square) value of j0 = I

A.

RMS(j0) =1√2· 100 A

4 cm2≡ 1.77 · 105 A

m2(84)

(a) stranded winding (b) massive winding

Figure 14: Detailed Current Density jNote: Due to the prescribed value the distribution of current and voltage excitation is equal.

38

5 Validation

In Figure 14b the distribution of a massive winding is visualized, which changes alongan arbitrary cross-section. The inner regions show a higher current density, whereasthe total integrated current leads to the same value of 100 A.

Although all models have different excitations, the current density leads to good ac-cordance in all four windings. Especially every current excitation can be replaced by aprescribed, precomputed voltage so that the same current density results. Further, thecase shows that the use of multiple windings in arbitrary positions inside of one modelcan be performed. Therefore the results in all windings are expectable and correct ina physical sense.

5.2.2 Compare Circular and Spiral Winding Model

In this comparison, the implemented winding model is compared to the already existingsolver HADAPT. Therefore two similar models with a bar, a winding, and an additionalwall are observed in order to verify the winding approach.

(a) circular winding (b) spiral winding

Figure 15: Geometry of Circular and Spiral Model

Figure 15 shows the two geometries, whereas the left model contains a full circularwinding and is simulated with the implemented winding model. The right modeldepicts a similar case, with a curved, spiral bar instead of the closed winding, which isexecuted with the already existing version. The winding has an inner radius of 5 cmand the wall a surface of 8 cm × 8 cm. The distance between the wall and the windingis 2 cm, and all components have a thickness of 1 cm. Both models contain a current

39

5 Validation

excitation of 100 A and 0 phase shift in the bar and the winding. Although thereare small differences in the model geometry the results, especially on the wall due toinduction effects, can be compared.

The material is set to copper for all conductors in both models with the parameterschosen in Section 5.2.1. The simulation is executed on a fine mesh with about 250000elements and with two frequencies, 50 Hz, and 0 Hz. For a better comparison, theinduced results of the current density j are observed on the wall for both almostidentical geometries. The results for the frequency 50 Hz are depicted in Figure 16.

(a) full winding (b) pseudo winding

Figure 16: Current Density j on Wall for 50 Hz

The visualization shows the induced current density on the wall, which is almost iden-tical for both models and confirms the embedding of the winding model in the existingsolver. Due to the excitation in the bar, the distribution is not exact symmetric to thehorizontal. For 0 Hz the induction effects are not present so that the resulting cur-rent density on the wall end up with zero in both geometries. Altogether the internalcomparison shows expectable results in all executed simulations.

5.3 External Comparison With Commercial Solver

In this Section, the winding model is compared to an external electromagnetic solvercalled MagNet. This solver provides a similar winding model with a simplification ofhollow cylinders. Therefore two different comparisons are performed in the following.

5.3.1 Massive Winding

In this validation, the commercial solver is used to show the accordance for the windingmodel, especially for massive windings. Instead of neglecting all induction effects as in

40

5 Validation

the comparison in Section 5.2.2, the skin effect is discussed with the aid of a massivewinding with one turn and an applied current excitation of 100 A. The turn has thesame dimensions and material parameters as the depicted turns in Section 5.2.1

The intensity of the skin effect depends on the frequency because for higher frequenciesthe skin effect increases and the current density is therefore concentrated at the edges.On the other hand, a frequency of 0 Hz leads to no induction effects. Hence the currentdensity is observed in a tangential component of an arbitrary cut surface, whereby thedistribution of the real part is depicted in Figure 17 for different frequencies.

(a) 50 Hz (b) 500 Hz (c) 5000 Hz

Figure 17: Current Density Re(j)z on Cross-Section for Different FrequenciesNote: The inner wall of the winding is the right edge of the depicted cross-section. The normal

component to the plane is the z-component.

In the case of 50 Hz, the induction effects are small, and the current density resultsin almost the prescribed distribution with a higher current at the center. For higherfrequencies, the skin effect is noticeable, so that the flow concentrates at the outerregions. For an extremely high voltage of 5000 Hz, the current is fully accumulatedat the right corners, which are the inner edges of the entire winding. Because of theprescribed current progression, the total current distribution is always unsymmetricand higher on the right side, which is a typical distribution for massive windings.

In order to verify the winding model, a comparison to the solver MagNet is performedfor all shown frequencies. The set mesh control for HADAPT results in a finer dis-cretization than the used MagNet mesh. The different discretizations and further theresults of both solvers for the 500 Hz case can be observed in Figure 18. For a bet-ter comparison, the results of the real- and imaginary part are depicted again for thez-component.

The shown results are almost identical in magnitude and distribution for both solversso that the winding model for massive windings shows excellent accordance with thesolver MagNet.

41

5 Validation

(a) real part HADAPT (b) real part MagNet

(c) imaginary part HADAPT (d) imaginary part MagNet

Figure 18: Current Density on Cross-Section for 500 Hz and Both SolversNote: The inner wall of the winding is the right edge of the depicted cross-section. The normal

component to the plane is the z-component.

5.3.2 Realistic Transformer

In the last validation, a realistic transformer is modeled and solved with both solvers,whereas the geometry is depicted in Figure 19a. The transformer contains three phaseswith outer primary and inner secondary coils, which end up with six coils in the inside.These coils - with up to thousand turns - are now modeled with the new implementedwinding model as hollow cylinders.

The prescribed excitation of the primary windings is 10 kA, whereas the phase angleshift varies for the different phases. For a three-phase current, the phase shift is 120,which leads to the different excitations in Figure 19b. Each primary winding consistsof 1000 turns so that with the amount of secondary turns the current Isecondary in the

secondary winding can be defined as Isecondary =Nprimary

Nsecondary· Iprimary = 50 kA. The phase

42

5 Validation

(a) geometry

Material Parameters

Part σ µ ε

Winding 5.96e+7 1.0 1.0Core 0.0 1000.0 1.0

Housing 5.0e+5 1.0 1.0

Note: σ in [ Sm ] and µ, ε in [−]

Excitation Values

Part Current Phase

Primary 1 10000 A 0

Primary 2 10000 A 120

Primary 3 10000 A 240

Secondary 1 50000 A −90

Secondary 2 50000 A 30

Secondary 3 50000 A 150

(b) materials and excitations

Figure 19: Simulation Setup for Transformer Model

angle shift between primary and secondary windings is constant -90. The resultingcurrent density is visualized for a frequency of 60 Hz and all windings in Figure 20.

(a) HADAPT

(b) MagNet

Figure 20: Current Density Re(j)x on Windings for HADAPT and MagNetNote: The visualization of the x-component allows a conclusion of the current direction.

43

5 Validation

The defined excitation and especially the phase shift of all secondary windings leadsto opposite directions of the current in the secondary windings. As seen in Figure 20the solution of both solvers are expectable similar in magnitude and distribution.

In addition to the six excited windings, a linear magnetic steel core exists beside thewindings. Due to an electric conductivity of zero the electric properties are neglectedinside of the core. Furthermore, an additional housing embraces the entire transformer.The entire geometry has the outer dimensions of 2 m × 2 m × 3 m and is depicted inFigure 19a. The material parameters for the different model parts are summarized inFigure 19b, where the windings and the housing are considered as copper and the coreas linear magnetic steel.

Obviously, the prescribed excitation induced a current and therefore losses on thetransformer housing, which is a typically observed phenomenon in business simulations.Due to the different excitations and positions of the windings, the resulting currentdensity on the housing is not homogenous. Figure 21 shows the distribution of theinduced current density j on the copper housing for both solvers.

(a) HADAPT (b) MagNet

Figure 21: Current Density j on Housing for HADAPT and MagNetNote: RMS values - root mean square values

The visualization of the current density shows that the distribution varies and is nothomogeneous. But in general, the comparison of both observed solvers leads to goodaccordance and almost the same results. With this shown property the losses on thehousing and the other transformer parts are observed and depicted in Table 2.

The resulting losses of all parts are due to the same excitations almost identical for both

44

5 Validation

Ohmic Losses Core Housing Primary Coils Secondary Coils

HADAPT 0 W 2.19e+9 W 5.271e+7 W 3.765e+7 WMagNet 0 W 2.23e+9 W 5.269e+7 W 3.764e+7 W

Table 2: Induced Losses in Transformer Model

electromagnetic solvers. Even the observed losses on the housing are acceptable witha deviation of 1.7 %. Because of a nonconductive core in the inside of the transformerand therefore no current, the losses are zero in this component.

Instead of the current density the magnetic field B on the non conductive steel core isobserved in a last comparison. The results are shown in the visualization in Figure 22.

(a) HADAPT (b) MagNet

Figure 22: Magnetic Field B on Core for HADAPT and MagNet

The comparison shows that the distribution and also the magnitude of the magneticfield lead to good accordance of both solvers. The maximum is reached at the edgesnext to the windings, whereas the minimum exists at the outer corners of the core.

Therefore all comparisons show accurate and expectable results so that the imple-mented winding model is a reasonable simplification for coils with a vast amount ofturns.

45

6 Conclusion

6 Conclusion

Final RemarksThe implemented winding model simplifies the simulation of coils by approximatingevery single turn as a coherent, hollow cylinder, whereby this approximation is applica-ble for an uniform, dense turn distribution. The prescribed excitation value is adaptedto the cylinder, so that almost the same results as in a circular coil are predictable.

The prescribed excitation value leads to an additional current density in the winding.Beside voltage or current excitation, there is a further distinction between strandedand massive windings, leading to different equation systems. Relevant in the caseof stranded winding is that the used winding tangential is divergence-free, which isprovided by solving a previous, smaller Poisson problem for all winding nodes. Withthis calculated winding tangential the full Maxwell equation system can be solvedwith a suitable solver. Of course, also the usage of both, voltage and current excitedwindings, inside of one model is possible.

For massive windings, every turn is modeled as a solid, where induction effects are nolonger neglected and therefore an inhomogeneous current distribution results. Insteadof solving a previous equation system a scaling of the resulting solving potential isnecessary, which is embedded inside of the equation system.

Beside some relative error calculations to the analytical solution and the simulation ofsimple winding examples, a simulation with a real industrial transformer is observedand compared to an external simulation tool. All validations for different configurationslead to good accordance so that the new model can be used in the following as a partof the electromagnetic solver HADAPT.

Future WorkThe new implemented winding model increases the productivity of the solver andenables the convenient simulation of complex transformer simulations, with arbitrarynumber of turns. For the continuous extension of the solver, some additional featuresare plannend to be realized inside of the solver.

Multiple windings with current excitation inside of one model lead to mutual influenceof the induction voltage. This influence can be calculated with a symmetric inductancematrix and is given by Ui = −

∑Nn=1 Lij ·iωIj. Row i of matrix Lij can be computed by

setting Ii = 1 and all other current zero. This results with the relation Lij =−Uj

iωIjto the

coefficients of the induction matrix. The realization of this requires the LES-solvingN times and is intended to come next. [8]

Due to acting forces, a mechanic coupling algorithm has to be applied, which calculates

46

6 Conclusion

the deformations of the modeled geometry [11]. Besides, new models for the treatmentof magnetic hysteresis are proposed. Further, also materials with anisotropic electricconductivities should be considered. These features require new approaches, whichhave to be integrated into the existing solver structure. [12]

At the moment two different coupling methods are provided for electrothermal compu-tations. The build-in thermal-mode calculates the temperature distribution with thepure heat equation in a postprocessing step, whereas the second coupling possibilityuses the external solver FLUENT. To improve the usability of the coupling the license-free tool OpenFOAM [9] can be used for an integrated two-way coupling. In such acoupling the thermal influences on the electric conductivity are considered so that atwo-step coupling loop can be established.

Altogether the solver is already applicable and used for real business simulations. Thenew implemented winding model provides the possibility of modeling entire transform-ers in reasonable computation times and offers a convenient interface to a CAD- andvisualization tool. Hence the presented winding model inside of the solver shows acomfortable method for simulating electromagnetic phenomena with windings.

47

Glossary

Glossary

B Magnetic flux [T]

H Magnetic field [A m−1]

D Electric displacement field [F V2 m−1]

E Electric Field [V m−1]

µ Relative permeability

σ Electric conductivity [S m−1]

ε Relative permittivity

R Resistance [Ω]

ρ Charge density [C m−3]

j Current density [A2 m−1]

ji Induced current density [A2 m−1]

j0 Applied current density [A2 m−1]

A Vector potential

ϕ Scalar potential

ϕi Scalar induced potential

ϕ Component of the induced potential ϕi

ψ Component of the induced potential ϕi

ϕ0 Scalar applied potential

Vring Winding Volume [m3]

Ω Winding domain in R3

∂Ω Surface of domain Ω

A Arbitrary cross-section of a winding [m]

Γ1 Surface on the left domain of the winding area

Γ2 Surface on the right domain of the winding area

48

Glossary

ΓD Subspace of Γ2 with defined Dirichlet condition

ω Angular velocity

Φ Prescribed phase angle

I Current [A]

Ii Induced current [A]

I0 Applied winding current [A]

I0 Prescribed peak value of the applied current I0 [A]

U Voltage [V]

Ui Induced voltage [V]

U0 Applied winding voltage [V]

U Prescribed peak value of the voltage U [V]

eϕ Winding tangential

eϕ Unit winding tangential with magnitude of one

~xwinding Midpoint of the winding

~xarea Midpoint of the cross-section

~rref Reference radius to midpoint of area)

~r Radius of a point to the winding midpoint

~rproj Projected radius of a point in the winding plane

r Scalar radius of a point, assumed as the projected radius

~a Axis of winding

~n Normal vector to the cross-section

~p Arbitrary point in the winding

Θ Angle in the winding plane

H(curl,Ω) v ∈ L2(Ω) ; curl(v) ∈ L2(Ω)

H(U) ψ ∈ H1(Ω) : ψ = 0 on all outer airbox surfaces

H1e (Ω) v ∈ H1(Ω) : v = const on all solid surfaces ;v = 0 on all airbox surfaces

49

Glossary

H1(Ω) v ∈ L2(Ω) ; Dα(v) ∈ L2(Ω) ∀ |α| ≤ 1

V v ∈ H(curl,Ω) : curl(vt) = 0 on ∂Ω,∫τv dS = 0

vt Tangential component of v

Vi Vectorial shape function

Wi Nodal shape function

Φi Nodal shape function for winding tangential system

ϕ′ Test function of potential ϕ

ψ′ Test function of potential ψ

A′ Test function of potential A

Nqp Number of quadrature points

ωk Quadrature weights at quadrature point

F Force acting on a winding [N]

c Scaling factor for massive windings

Ascaled Scaled vector potential A

ϕi, scaled Scaled potential ϕi

ϕ0, scaled Scaled potential ϕ0

[·] Jump of the potential ·

Iscaled Scaled total current [A]

Iwanted Prescribed current [A]

δϕ0 Relative error in L2-norm of ϕ0

δeϕ Relative error in L2-norm of eϕ

50

References

References

[1] F. Assous, P. Degond, E. Heintze, P.A. Raviart, and J. Segre. On a Finite-Element Method for Solving the Three-Dimensional Maxwell Equations. Journalof Computational Physics, pages 223–237, 1993.

[2] Raffael Casagrande. Implementation of impedance boundary conditions. 2012.

[3] Richard W. Cottle. Manifestations of the Schur complement. Rendiconti delSeminario Matematico e Fisico di Milano, pages 189–211, 1975.

[4] Wolfgang Dahmen and Arnold Reusken. Numerik fur Ingenieure und Naturwis-senschaftler. 2008.

[5] Thomas Friedrich Denny. Losen grosser dunnbesetzter Gleichungssysteme uberendlichen Primkorpern. PhD thesis, Universitat des Saarlandes, 1997.

[6] Ralf Hiptmair. Finite elements in computational electromagnetism. Acta Numer-ica, pages 237–339, 2002.

[7] Ralf Hiptmair, Florian Kramer, and Jorg Ostrowski. A robust Maxwell formula-tion for all frequencies. IEEE Transactions on Magnetics, pages 682–685, 2008.

[8] Ralf Hiptmair, Jorg Ostrowski, and Florian Kramer. HADAPT Full Maxwell(RMF-module). Technical report, 2017.

[9] Hrvoje Jasak, Aleksandar Jemcov, and Zeljko Tukovic. OpenFOAM : A C ++Library for Complex Physics Simulations. International Workshop on CoupledMethods in Numerical Dynamics, pages 1–20, 2007.

[10] Benjamin Menkuc and Ralf Wiechmann. Herleitung der Divergenz in Zylinderko-ordinaten ausgehend von kartesischen Koordinaten. Technical report, 2004.

[11] Philipp Muller. Macroscopic Electro Thermal Simulation of Contact Resistances.2016.

51

References

[12] PLCRC Research Center ABB. Simulation Toolbox Website. URL http:

//simulation-toolbox.abb.com. (retrieved on 2018-02-23).

[13] H.T. Rathod, K.V. Nagaraja, B. Venkatesudu, and N.L. Ramesh. Gauss Legendrequadrature over a triangle. J. Indian Inst. Sci. 84, pages 183–188, 2004.

[14] Peter Russer, Mauro Mongiardo, and Leopold B. Felsen. Electromagnetic FieldComputation by Network Methods. 2009.

[15] Youcef Saad. SPARSKIT: A basic tool kit for sparse matrix computations. Tech-nical Report NASA, 1990.

[16] Olaf Schenk, Klaus Gartner, and Wolfgang Fichtner. PARDISO: a high-performance serial and parallel sparse linear solver in semiconductor device simu-lation. Future Generation Computer Systems, pages 69–78, 2001.

[17] C. F. Smith, A. F. Peterson, and R Mittra. The biconjugate gradient methodfor electromagnetic scattering. IEEE Transactions on Antennas and Propagation,1990.

[18] A. Taflove and M.E. Brodwin. Numerical Solution of Steady-State Electromag-netic Scattering Problems Using the Time-Dependent Maxwell’s Equations. IEEETransactions on Microwave Theory and Techniques, pages 623–630, 1975.

[19] H. A. van der Vorst. Bi-CGSTAB: A Fast and Smoothly Converging Variant ofBi-CG for the Solution of Nonsymmetric Linear Systems. Society for Industrialand Applied Mathematics, 1991.

[20] Thomas Weiland. Die Diskretisierung der Maxwell-Gleichungen. Phys. UI. 42(1986) Nr. 7, pages 191–201, 1986.

[21] H.A. Wheeler. Formulas for the Skin Effect. IEEE RFIC Virtual Journal IEEERFID, 1942.

52

http://simulation-toolbox.abb.com

http://simulation-toolbox.abb.com

FEM Implementierung von Spulen Modellen f ur Transformator ... · Transformator Simulationen FEM...

Documents

Transcript of FEM Implementierung von Spulen Modellen f ur Transformator ... · Transformator Simulationen FEM...