Herbert De

Ruth Sabariego

Mr Ronan

I take advantage in the following of these different paragraphs to thank all the people who

Introduction

The 19th century saw the birth of the theory of electromagnetism with the great pioneers of this whole branch of modern physics. Around 1820, Oersted and Ampère were the first to demonstrate the link between magnetism and electricity, which had long been studied as two separate phenomena. Subsequent work by Faraday and others led to Maxwell's dynamic theory of the electromagnetic fields, published in 1864. This theory was validated in 1888 by Hertz with the discovery of electromagnetic waves propagating at the speed of light.

However, depending on the frequency of the studied phenomena, different problems can be considered. The first one is focused on the high frequency cases such as hyperbolic problems. The second one is applied to the low frequency cases, in particular when the inductive and capacitive phenomena are decoupled. Nevertheless, the limit between high and low frequencies is not well defined. Indeed, for certain socalled intermediate frequencies, the radiation effects can be neglected. Besides, the coupled capacitive-inductive effects should be taken into account since the resonance phenomena may appear.

Neglecting the radiation effects, several approximated models can be derived from Maxwell's equations known as static or quasistatic models. In the literature, the electrostatic, electrokinetic, and magnetostatic problems are widely used in low frequency to capture the decoupled capacitive, resistive, and inductive effects, respectively [Haus and Melcher, 1989a;[START_REF] Badics | Full wave potential formulation with low-frequency stability including ohmic losses[END_REF].

Moreover, when the frequency of the voltage or the current supply applied to windings of electrical devices increases, two different phenomena cannot be neglected.

On the one hand, the electric field in the regions close to the windings becomes im-INTRODUCTION portant, in this case, one special attention should be taken on the coupled resistivecapacitive effects. The electro-quasistatic (EQS) model [START_REF] Steinmetz | Electro-quasistatic field simulations based on a discrete electromagnetism formulation[END_REF] is commonly used to fields resulting from high voltage, withstand applications or microelectronics. On the other hand, the skin effect appears on the distribution of the current density in the winding, the coupled resistive-inductive effects have to be taken into account where the magneto-quasistatic (MQS) [Haus and Melcher, 1989a] model is widely used.

Besides, if the two phenomena should be characterized at the same time, the full

Maxwell model [START_REF] Hiptmair | A robust Maxwell formulation for all frequencies[END_REF][START_REF] Badics | Full wave potential formulation with low-frequency stability including ohmic losses[END_REF] should be used even it is more time-consuming and has an instability issue in low frequency as reported in [START_REF] Hiptmair | A robust Maxwell formulation for all frequencies[END_REF][START_REF] Bebendorf | Hierarchical matrix preconditioning for lowfrequency-full-Maxwell simulations[END_REF]. If the operation frequency is considered in the intermediate range, where the coupled capacitive-inductive effect can be handled while the radiation effects are neglected, the Darwin model can be used as a surrogate model of the full Maxwell system [START_REF] Larsson | Electromagnetics from a quasistatic perspective[END_REF][START_REF] Fang | Darwin approximation to Maxwell's equations[END_REF][START_REF] Koch | A low-frequency approximation to the Maxwell equations simultaneously considering inductive and capacitive phenomena[END_REF][START_REF] Eller | A symmetric low-frequency stable broadband Maxwell formulation for industrial applications[END_REF][START_REF] Zhao | A novel gauged potential formulation for 3-D electromagnetic field analysis including both inductive and capacitive effects[END_REF][START_REF] Clemens | A Darwin time domain scheme for the simulation of transient quasistatic electromagnetic fields including resistive, capacitive and inductive effects[END_REF].

The Darwin approximation is proposed around 1920 by Charles Galton Darwin, who is the grandson of Charles Darwin. It has been proved in 1992 by [START_REF] Degond | An analysis of the darwin model of approximation to maxwell's equations[END_REF].

The main objective of the thesis is to implement tools in the code_Carmel (Code Avancé de Recherche en Modélisation Électromagnétique) in order to quantify the electric field in dielectrics, in particular the Darwin model. On the other hand, before modeling electrical machines, the Darwin model should be studied as well as its limits. In the literature, the classification or the theoretical limitations of the different quasistatic models are mainly based on the theoretical considerations, in particular, with several assumptions for the computational domain. For example, only the conductive domain is considered, not for the case with multi-domain. The existing reference [START_REF] Koch | Different types of quasistationary formulations for time domain simulations[END_REF] showed a comparison of the results obtained with different quasistatic models only in the time domain while a simple axisymmetric test model represented by a parallel plate capacitor is used. The comparison results were illustrated for the electric field. Our motivation here is to investigate the limits of the different quasistatic models in the general case [START_REF] Taha | Numerical simulation-based investigation of the limits of different quasistatic models[END_REF], in particular, for the computational domain including both conductive and non-conductive domains, which has not been well addressed in the literature. To achieve this objective, in this work, we first validate the numerical results obtained by the Darwin model with measurements. Furthermore, to the best of the author's knowledge, the comparison of the Darwin model with real measurement results has not been reported in the literature.

In addition, we are interested in integrating special elements such as shell elements [START_REF] Ren | Degenerated Whitney prism elements-general nodal and edge shell elements for field computation in thin structures[END_REF] in the Darwin model in order to further reduce the computational cost. In this thesis, we develop and implement shell elements for accurately considering thin-conducting pads, traces, layers in printed circuit boards (PCBs) while avoiding their volume mesh. This proposed approach is successfully validated on two industrial PCBs.

The manuscript of this thesis is divided into four chapters. In the first chapter, which is entitled "Modeling of electromagnetic problems", a brief introduction of different mathematical models resulting from the Maxwell equations is given, with detailed assumptions. The associated numerical method as the finite element method is recalled, with the introduction of the continuous variational formulations, discrete variational formulations. Then, the corresponding formulations for each quasistatic models, namely EQS, MQS, and Darwin are given.

The

Maxwell's equations

More and more electromagnetic devices are subjected to high frequency excitations.

These operating conditions contribute to the development of undesirable effects in magnetic cores and wound parts (capacitive effects between turns, conductive losses, etc.). The modeling of high frequency devices requires a complete solving of Maxwell's equations, taking into account in particular the displacement currents.

The modeling of this magnetic-electric interaction should then allow a better understanding of the real functioning of these electromagnetic devices. Maxwell's equations constitute a system of partial differential equations, they describe the electromagnetic behavior of any system. Let Ω ⊂ R 3 be a domain of study with Γ the boundary. Ω c ⊂ Ω is the sub-domain for conductive part and Ω nc is a subdomain for non-conductive part. In continuous domain, electromagnetic phenomena are described by four functions which depend on time t ∈ [0, T ] where T denotes a positive real, and space coordinates x, with values in R 3 :

curl H(x, t) = J(x, t) + ∂D ∂t (x, t), Maxwell -Ampere (1.1) curl E(x, t) = - ∂B ∂t (x, t), Maxwell -Faraday (1.2) div B(x, t) = 0, Maxwell -Thomson (1.3) div D(x, t) = ρ(x, t), Maxwell -Gauss (1.4)
where

• H is the magnetic field [A.m -1 ],
• J is the current density [A.m -2 ],

• D is the electric flux density [C.m -2 ],

• E is the electric field [V.m -1 ],

• B is the magnetic flux density [T],

• ρ is the volume charge density [C.m -3 ].

The conduction current density J can represent two terms: the induced current density J ind in conductive materials and J s the uniform current density in multifilament inductors.

ELECTROMAGNETIC MODELING

In the following, we write Maxwell's system directly as

curl H = J + ∂D ∂t , (1.5) curl E = - ∂B ∂t , (1.6) div B = 0, (1.7) div D = ρ.
(1.8)

The equations (1.5)-(1.8) have a purely geometric character and are independent of the characteristics of the domains. To model the behavior of the materials, it is necessary to add the constitutive laws or constitutive relations. They express the quantities D, B and J as functions of the electromagnetic fields E and H.

Behavior laws

The fields in Maxwell's equations are further related to each other by the material relations. These material equations are known as constitutive relations. When the material is insulating, σ tends to zero, therefore we have J = 0, which indicates that is no current flowing in the material. In the case of a perfect conductor, σ is infinity which equivalent to E = 0 and H = 0. In perfect material, that is to say the material for which the constitutive laws are linear, the following relations are verified:

Electric constitutive law

Let us denote ε the electric permittivity, where ε = ε r ε 0 with ε 0 the permittivity in vacuum equal to ε 0 = (36π10 9 ) -1 m -3 .kg -1 .s 4 .A 2 and ε r the relative permittivity.

The constitutive electric law depends on the electric field. Then, the dielectric constitutive relation between D and E is written as follows

D = εE = ε 0 ε r E.
(1.9)

Magnetic constitutive law

Let us denote µ the magnetic permeability, where µ = µ r µ 0 with µ 0 the permeability in vacuum equal to 4π10 -7 m.kg.s -2 .A -2 and µ r the relative permeability. The linear magnetic constitutive law [START_REF] Jackson | Classical Electrodynamics[END_REF] is written as follows

B = µH = µ 0 µ r H.
(1.10)

The equation (1.10) describes the behavior of materials for the applications targeted in electrical engineering. For ferromagnetic materials, relatively complex models can CHAPTER 1. MODELING OF ELECTROMAGNETIC PROBLEMS be used which take into account the phenomenon of hysteresis [START_REF] Jiles | Theory of ferromagnetic hysteresis[END_REF]].

In the case of hard materials (permanent magnets), the phenomenon of remanence is introduced, and we obtain a law of the form:

B = µ 0 H + B r . (1.11)
where B r is the remanent induction.

Ohm's law

Similarly, by denoting σ expressed in [S.m -1 ] the electric conductivity, Ohm's law [START_REF] Jackson | Classical Electrodynamics[END_REF] comes to complement the previous two laws. This law was discovered by George Simon Ohm in 1827, allows to link the current density and the electric field. It is written as follows

J = σE.
(1.12)

The electric conductivity law (1.12) is often linear but in semi-conductors, this law can be non-linear. In addition, the conductivity depends strongly on the temperature. Certain simulations, in particular those intended to improve the efficiency of electrotechnic systems, need to take this effect into account. In this case, the temperature distribution in the conductive parts must be determined. Then, multiphysics modeling is necessary by coupling an electromagnetic model and a thermal model.

Boundary conditions

In addition to the Maxwell system equations (1.5)-(1.8) and the constitutive laws (1.9), (1.10), and (1.12), it is necessary to define the boundary conditions on the fields in order to obtain a unique solution to the problem. In this part, we present the homogeneous boundary conditions of magnetic and electric quantities. It is assumed in the following that the electric charge density ρ is zero.

Magnetic quantities

H × n | Γ H = 0, (1.13a) B • n | Γ B = 0 (1.13b) 1.1. ELECTROMAGNETIC MODELING Figure 1.1: Studied domain.
where n is a unit vector to Γ and oriented towards the outside of Ω.

The boundary condition defined in (1.13a) means that the tangential component of the magnetic field H is zero and that it is normal to the Γ H boundary. In the same way, the equation (1.13b) means that the normal component of the magnetic flux density B is zero, then B is tangent to the boundary Γ B .

Electric quantities

The conductive domain Ω c ⊂ Ω having a boundary which has a common part with Γ and another included strictly in Ω denoted Γ J . The boundary conditions of electric quantities are such that

E × n | Γ E = 0, (1.14a) J • n | Γ J = 0, (1.14b) D • n | Γ D = 0. (1.14c)
The boundary conditions (1.14b) and (1.14c) are obtained by combining (1.13a) with (1.5), while the condition (1.13b) is obtained by combining (1.14a) with (1.6).

Approximation of quasistatic regimes

The solving of the complete model of Maxwell's equations is not systematic, depending on the applications, certain terms are negligible on the scales of physical quantities.

Neglecting the radiation effects, several approximated models can be derived from Maxwell's equations based on the known Galilean limits [START_REF] Haus | Electromagnetic Fields and Energy[END_REF][START_REF] Rapetti | Implications of Galilean electromagnetism in numerical modeling[END_REF]. There are three main different approaches, one is

∂E ∂x ∼ E ′ L and ∂E ∂t ∼ E ′ T . (1.15)
In vacuum, by the Faraday's law equation (1.6), we have

E ′ L ∼ B ′ T → E ′ ∼ |v|B ′ . (1.16)
Similarly, the Ampere's law equation (1.5) is written as

B ′ µL ∼ εE ′ T → B ′ ∼ |v| c 2 E ′ .
(1.17)

We denote the notation "x ∼ y" means that the two quantities x and y have the same order. In the quasistatic regime or Galilean regime, the two equations namely Faraday's law and Ampere's law cannot be coupled. Suppose that these equations are coupled, the two ordering defined in (1.16) and (1.17) are different and valid on the same time, then if one replaces the other assumption, we get |v|∼ c, that is a contradiction with the Galilean regime assumption |v|≪ c. This proof verifies that these two orders are not validated simultaneously.

Suppose that the assumption (1.17) holds and E ′ ≫ |v|B ′ , the EQS approximation has to be adopted. Then, the term ∂B ∂t defined in Faraday's law as given in the equation (1.6) becomes negligible, which implies that the induced current density is

1.1. ELECTROMAGNETIC MODELING neglected. The EQS model is given as follows div (J + ∂D ∂t ) = 0, curl E = 0, div J = 0, div D = 0.
Suppose that the assumption (1.16) is hold and B ′ ≫ (|v|/c 2 )E ′ , the MQS model has to be adopted. Then, the term ∂D ∂t defined in the Ampere's law as given in the equation (1.5) becomes negligible, which implies that the displacement current density is neglected in the dielectrics. The MQS model is given as follows

curl E = - ∂B ∂t , curl H = J, div B = 0, div J = 0.
To obtain the Darwin model, based on the Helmholtz decomposition [START_REF] Jackson | Classical Electrodynamics[END_REF],

the electric field E is splitted into two parts, an irrotational part E irr , which is curl free and a solenoidal part E sol , which is divergence free

E = E irr + E sol , curl E irr = 0, div E sol = 0. (1.20)
We rewrite the Maxwell equations using this decomposition in vacuum

curl E sol = - ∂B ∂t , (1.21) curl H = ∂ ∂t (εE irr ) + ∂ ∂t (εE sol ), (1.22) div B = 0, (1.23) div (εE irr ) = 0. (1.24)
Based on the above equations, there exists two orders of magnitude [START_REF] Kruger | The three quasistatic limits of the Maxwell equations[END_REF]]

E ′ sol B ′ ∼ |v| and E ′ irr B ′ ∼ 1 |v| . (1.25)
curl H = J + ∂ ∂t (εE irr ), (1.27) div B = 0, (1.28) div D = 0, (1.29) div J = 0.
(1.30)

Different models

The electromagnetic models are derived from simplifying assumptions, while the main question is: what are the parameters that allow these assumptions to be applied? In reality, the physical phenomena presented in an electromagnetic system are of three types: resistive, inductive, and capacitive. Also, the choice is dictated by the presence or not as well as the preponderance over others of these phenomena in the device considered.

Choice of a model and associated physical phenomena

Table 1.1: Electromagnetic models classified according to the physical phenomena taken into account.

We assume that the materials are linear and homogeneous. Then, based on Maxwell's 1.2. DIFFERENT MODELS system (1.5)-(1.8), the Full-wave formulations are written

curl (curl E) + ∂(µσE) ∂t + ∂ 2 (µεE) ∂t 2 = 0, (1.31) curl (curl H) + ∂(µσH) ∂t + ∂ 2 (µεH) ∂t 2 = 0. (1.32)
To define the quantitative criteria making it possible to choose the formulations of this Table 1.1. We write the equations (1.31) and (1.32) in the frequency domain [START_REF] Geuzaine | Numerical techniques for forward modelling[END_REF][START_REF] Parent | Apport de la simulation numérique à l'amélioration de l'efficacité énergétique et du système d'isolation électrique des dispositifs électrotechniques[END_REF] curl

(curl x) + jωµσx -ω 2 µεx = 0 (1.33)
where x represents the electric or magnetic fields and ω its pulsation. Define δ the skin thickness, L the domain size, and λ the wavelength such as

δ = 2 σµω and λ = 2πc ω .
(1.34) Thus, the relation (1.33) reads as follows

3 L 2 - 2i δ 2 + 4π 2 λ 2 x = 0. (1.35)
In consideration of these relations, it is possible to define three parameters and choose our model according to them. We have

k 1 = λ L 2 , k 2 = δ L 2 and k 3 = λ δ 2 (1.36)
where the term k 1 allows us to specify if the propagation effects are dominant or negligible, k 2 allows to determine if the size of the system allows the current to flow with or without skin effect and k 3 allows to know which one of conduction or displacement currents is negligible, or equivalent, with respect to each other. Finally, the following classification can be concluded:

k 1 ≲ 1 : Full-wave;

-k 3 ≈ 1 : Electro-quasistatic; -k 2 ≲ 1 : Magneto-quasistatic;
k 1 ≫ 1 : Electric and magnetic phenomena are decoupled.

Quasistatic models on potentials formulation

In low frequencies, the inductive effects associated with eddy currents constitute the major effect on the functionality of the corresponding electrotechnical devices. However, especially for inverter powered electric machines, where the machine windings are exposed to higher harmonic field components from the fast switching inverter, the capacitive effects also become relevant. So to properly characterize the behavior of an electric machine, it should be necessary to include capacitive effects in addition to inductive and resistive effects. The common mode voltage caused by the inverter in combination with the dielectric insulation of the machine windings generates common mode currents at higher frequencies [START_REF] Maki-Ontto | Induction motor model for the analysis of capacitive and induced shaft voltages[END_REF].

Many researchers have adopted field formulations to solve Maxwell's equations [START_REF] Bossavit | Two variational formulations of the three-dimensional eddycurrent problem are discussed and compared[END_REF][START_REF] Ren | A new hybrid model using electric field formulation for 3-D eddy current problems[END_REF]Albanese and Rubinacci, 1990;[START_REF] Daveau | New (e,h) formulation coupling a finite element method and a boundary integral method for the computation of the interaction of waves with a conducting domain[END_REF][START_REF] Bossavit | Mixed-hybrid methods in magnetostatics: complementarity in one stroke[END_REF]. However, we prefer the formulations in scalar or vector potentials for the solving of these equations. Indeed, the potentials are mathematical tools introduced to facilitate the solving. Depending on the physical significance of the potentials, we distinguish two groups of formulations: electric formulations and magnetic formulations. We define φ and A the electric scalar potential and the magnetic vector potential, respectively.

To introduce the Full-wave potential formulation, A and φ are considered. Since Ω is a simply connected domain and according to (1.7), B is divergence-free, then there exists a magnetic vector potential (MVP) A such that B = curl A.

(1.37)

Likewise, the electric field E is expressed as a function of the potentials. We replace (1.37) in the equation (1.6), then

curl E + ∂A ∂t = 0. (1.38)
Therefore, there is an electric scalar potential φ defined in Ω such that

E = - ∂A ∂t -grad φ.
(1.39)

Due to the explicit appearance of the electric scalar potential (ESP) φ, the Maxwell Full-wave formulations using A and φ are very practical for handling voltage excitations by simply imposing Dirichlet conditions on the boundaries of conductor terminals. It should be mentioned here that the models to be studied in this thesis does not take into account the coils J s or the magnets B r .

DIFFERENT MODELS

Based on (1.5) and on the constitutive behavior laws, the potential formulation of the full Maxwell model reads as follows

curl 1 µ curl A + σ ∂A ∂t + grad φ + ∂ ∂t ε ∂A ∂t + grad φ = 0, (1.40a) div σ ∂A ∂t + grad φ + ∂ ∂t ε ∂A ∂t + grad φ = 0. (1.40b)
In order to model the complete behavior of an electromagnetic device, such as in the complicated complex conductive windings of high frequency transformers and electric machines, as well as fast transient overvoltages in traditional transformers due to lightning [START_REF] Ostrowski | Transient Full Maxwell Computation of Slow Processes[END_REF], it should be taken into account the capacitiveinductive coupled effects. Then, it is essential to integrate the displacement current in the numerical modeling.

In the following, the three different models namely EQS, MQS, and Darwin models issue from the full Maxwell system will be studied in terms of potentials.

Electro-quasistatic model

Electrodynamics is a means of studying electric effects in an unsteady state which means that the phenomena considered depend on time. The objective of the electrodynamic or electro-quasistatic (EQS) model [START_REF] Steinmetz | Electro-quasistatic field simulations based on a discrete electromagnetism formulation[END_REF]] is to study a device which takes into consideration two important physical phenomena: conduction currents and displacement currents. The similarity of the two quantities was found by Maxwell by supposing that they are both sources of a magnetic field H. Ω nc can be added to represent the dielectric part. Whereas, the EQS approximation of Maxwell's equations neglects inductive phenomena, while capacitive and resistive phenomena are taken into account. Electro-quasistatic fields must be taken into account with capacitive-resistive structures such as, for example, resistance coated insulating structures and high voltage bushings or surge arresters used in high voltage electric technology. This model is valid when the characteristic length of magnetic phenomena is low w.r.t. the wavelength [START_REF] Steinmetz | Electro-quasistatic field simulations based on a discrete electromagnetism formulation[END_REF][START_REF] Weida | Design of zno microvaristor material stress-cone for cable accessories[END_REF]. Finite element formulations using standard numerical linear algebra schemes for this type of problem have already been proposed in [START_REF] Kurita | Dc flashover voltage characteristics and their calculation method for oil-immersed insulation systems in hvdc transformers[END_REF][START_REF] Preis | Numerical computation of transient quasistatic electric fields[END_REF][START_REF] Preis | Fem simulation of thermistors including dielectric effects[END_REF].

We take the EQS system defined in sub-section 1.1.4 by adding the boundary conditions. Then, the EQS problem is given as follows div (J + ∂D ∂t ) = 0, (1.41a

)

curl E = 0, (1.41b) div J = 0, (1.41c) div D = 0, (1.41d) E × n | Γ E = 0, (1.41e) J • n | Γ J = 0, (1.41f) D • n | Γ D = 0. (1.41g)
φ formulation

The electric scalar potential φ can be derived from the Maxwell equation (1.41b)

such as E = -grad φ.

(1.42)

H × n | Γ H = 0, (1.48) E × n | Γ E = 0, (1.49) J • n | Γ J = 0. (1.50)
In general, we also say eddy current problem for the MQS problem. Magnetodynam- CHAPTER 1. MODELING OF ELECTROMAGNETIC PROBLEMS

A -φ formulation

We note here that the potential A is not unique and there exist an infinite values of magnetic vector potentials since A is defined up to a gradient. Several solutions exist to ensure the uniqueness of A, we cite here the Coulomb gauge (div A = 0) [START_REF] Bellac | Galilean electromagnetism[END_REF] and the tree gauge (A • w = 0) [Albanese and Rubinacci, 1990], where w represents an arbitrary vector field whose field lines connect all the points of the domain and do not close. This gauge consists in restricting the number of unknowns of A [START_REF] Carpenter | Comparison of alternative formulations of 3-dimensional magnetic-field and eddy-current problems at power frequencies[END_REF].

Follow the same procedure made for the Full-wave system, where the magnetic flux density B is written as a function of A (1.37) and the electric field E is written as a function of A and φ (1.39), the equations defined in (1.44)-(1.50) are written now as

curl 1 µ curl A + σ ∂A ∂t + grad φ = 0 on Ω, (1.51a) div σ ∂A ∂t + grad φ = 0 on Ω c , (1.51b) curl A × n = 0 on Γ H , (1.51c) A × n = 0 on Γ B (Γ E ), (1.51d
)

grad φ • n = 0 on Γ H , (1.51e) φ = k on Γ E . (1.51f)
The aboved mentioned system is called the MQS model in A -φ potential formulation. We should mention here that φ is defined only in the conductive sub-domain Ω c . This model does not take into account the capacitive effect when the frequency increases. In the following, the Darwin model which able to capture the coupled capacitive-inductive effects will be presented.

Darwin model

Darwin's approximation of Maxwell's equations is introduced to remove the most rigid time scale in electromagnetic simulations, which is called wave propagation, but it keeps an important part of physics, especially the phenomena at low frequencies.

In

H × n | Γ H = 0, E sol × n | Γ E = 0, J • n | Γ J = 0.

A -φ formulation

Similar to the Full-wave and MQS model, B and E can be expressed on potential functions where

B = curl A, E = - ∂A ∂t -grad φ.
From the Helmholtz decomposition (1.20), E irr can be represented by the gradient of the scalar electric φ while E sol can be represented by the time derivative of the magnetic vector

E = E irr + E sol = -grad φ - ∂A ∂t .
(1.58)

We replace now the expressions (1.37) and (1.39) in (1.53), we obtain

curl 1 µ curl A + σ ∂A ∂t + grad φ + ∂ ∂t (εgrad φ) = 0. (1.59)
The charge conservation equation is written in Ω with ρ = 0 as follows

div σ ∂A ∂t + grad φ + ∂ ∂t (εgrad φ) = 0. (1.60) CHAPTER 1. MODELING OF ELECTROMAGNETIC PROBLEMS
In conclusion, the A -φ formulation for the Darwin model reads

curl 1 µ curl A + σ ∂A ∂t + grad φ + ∂ ∂t (εgrad φ) = 0 on Ω, (1.61a) div σ ∂A ∂t + grad φ + ∂ ∂t (εgrad φ) = 0 on Ω, (1.61b) curl A × n = 0 on Γ H , (1.61c) A × n = 0 on Γ E , (1.61d) grad φ • n = 0 on Γ H , (1.61e) φ = k on Γ E . (1.61f)
There are ungauged formulations proposed in the frequency domain [START_REF] Koch | A low-frequency approximation to the Maxwell equations simultaneously considering inductive and capacitive phenomena[END_REF] and lately reported for the time domain [START_REF] Ho | Application of edge elements to 3-D electromagnetic field analysis accounting for both inductive and capacitive effects[END_REF] neglecting the second order time terms from (1.40) for low frequency applications, i.e. when the electric size of the device is much shorter than the smallest wavelength, and taking into account the coupled inductive and capacitive effects.

In the next part, the finite element method will be presented in order to solve the potential formulations defined in the sub-section 1.2.2.

Finite element method 1.3.1 Continuous modeling 1.3.1.1 Continuous functional spaces

The electromagnetic fields presented in Maxwell's equations are defined in mathematical spaces called "Sobolev spaces". The latter are vector spaces well suited to partial differential problems. Moreover, the use of the finite element method requires the introduction of certain functional spaces which one defines briefly in the following. Typically, spaces denoted by bold uppercase letters are vector field spaces, and spaces denoted by normal uppercase letters are scalar field spaces.

Lebesgue spaces

We start by introducing the "Lebesgue spaces", we denote by L 2 (Ω), the space of square scalar fields integrable on Ω and L 2 (Ω), the space of vector fields whose square of the Euclidean norm (denoted

|•|) is integrable on Ω, such as L 2 (Ω) = {u : Ω → R | Ω |u(x)| 2 dΩ < +∞}, (1.62) L 2 (Ω) = {u : Ω → R 3 | Ω |u(x)| 2 dΩ < +∞}.
(1.63)

1.3. FINITE ELEMENT METHOD
The respective dot products of L 2 (Ω) and L 2 (Ω) spaces are defined by

(u, v) L 2 (Ω) = Ω uv dΩ, ∀u, v ∈ L 2 (Ω), (1.64) (u, v) L 2 (Ω) = Ω u • v dΩ, ∀u, v ∈ L 2 (Ω) (1.65)
associated with the following norm

||u|| 2 L 2 (Ω) = (u, u) 2 L 2 (Ω) , ∀u ∈ L 2 (Ω), (1.66) ||u|| 2 L 2 (Ω) = (u, u) 2 L 2 (Ω) , ∀u ∈ L 2 (Ω).
(1.67)

Sobolev spaces

Similarly, the Sobolev spaces [START_REF] Monk | Finite Element Methods for Maxwell's Equations[END_REF][START_REF] Girault | Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms[END_REF]] H 1 (Ω) and

H 1 (Ω) are defined as H 1 (Ω) = {u ∈ L 2 (Ω) | ∂u ∂x i ∈ L 2 (Ω), ∀i = 1, • • • , m}, (1.68) H 1 (Ω) = {u ∈ L 2 (Ω) | grad u ∈ L 2 (Ω)}.
(1.69)

Functional spaces of differential operators

The functional spaces of the differential operators, the gradient, the rotational, and the divergence are defined in this section.

Partial differential equations are solved in the domain Ω and they involve differential operators such as gradient, rotational, and divergence. The spaces of these differential operators are defined such that grad, curl, and div belong to L 2 (Ω) or L 2 (Ω), such as [Bossavit, 1988a,b]:

H(grad, Ω) = {u ∈ L 2 (Ω) | grad u ∈ L 2 (Ω)} = H 1 (Ω), (1.70) H(curl, Ω) = {u ∈ L 2 (Ω) | curl u ∈ L 2 (Ω)}, (1.71) H(div, Ω) = {u ∈ L 2 (Ω) | div u ∈ L 2 (Ω)}. (1.72)
The dot products of these spaces (1.70), (1.71), and (1.72) are defined respectively from the dot products (1.64) and (1.65) such as

(u, v) H 1 (Ω) = (u, v) L 2 (Ω) + (grad u, grad v) L 2 (Ω) , ∀u, v ∈ H 1 (Ω), (1.73) (u, v) H(curl, Ω) = (u, v) L 2 (Ω) + (curl u, curl v) L 2 (Ω) , ∀u, v ∈ H(curl, Ω), (1.74) (u, v) H(div, Ω) = (u, v) L 2 (Ω) + (div u, div v) L 2 (Ω) , ∀u ∈ H(div, Ω) (1.75)
associated with the following norm

||u|| 2 H 1 (Ω) = ||u|| 2 L 2 (Ω) +||grad u|| 2 L 2 (Ω) , ∀u ∈ H 1 (Ω), (1.76) ||u|| 2 H(curl, Ω) = ||u|| 2 L 2 (Ω) +||curl u|| 2 L 2 (Ω) , ∀u ∈ H(curl, Ω), (1.77) ||u|| 2 H(div, Ω) = ||u|| 2 L 2 (Ω) +||div u|| 2 L 2 (Ω) , ∀u, v ∈ H(div, Ω).
(1.78)

In the following, we define three functional sub-spaces of the differential operators.

These functional spaces also satisfy the boundary conditions and are defined as [START_REF] Bossavit | Whitney forms: a class of finite elements for three-dimensional computations in electromagnetism[END_REF][START_REF] Bossavit | Électromagnétisme, en vue de la modélisation[END_REF]:

H 1 Γ (Ω) = {u ∈ H(grad, Ω), u = cst on Γ}, (1.79) H Γ (curl, Ω) = {u ∈ H(curl, Ω), u × n = 0 on Γ}, (1.80) H Γ (div, Ω) = {u ∈ H(div, Ω), u • n = 0 on Γ}. (1.81)
Consequently, a link exists between the spaces defined above grad, curl, and div called the Rham complex [START_REF] Bossavit | Whitney forms: a class of finite elements for three-dimensional computations in electromagnetism[END_REF]]

H 1 (Ω) grad ===⇒ H(curl, Ω) curl = == ⇒ H(div, Ω) div ==⇒ L 2 (Ω).
(1.82)

In fact, we have

∀f ∈ H 1 (Ω) ⇒ grad f ∈ H(curl, Ω), (1.83) ∀f ∈ H(curl, Ω) ⇒ curl f ∈ H(div, Ω), (1.84) ∀f ∈ H(div, Ω) ⇒ div f ∈ L 2 (Ω).
(1.85)

For the rest, we mention here that E and H are defined in H(curl, Ω) and J and B in H(div, Ω).

Continuous variational formulations EQS problem in φ formulation

Let φ ′ be a test function defined on Ω. By applying the Green formula (A.2) to

(1.43), the variational formulation associated with the problem (1.43) is written

Find φ ∈ H 1 (Ω) such as (σgrad φ, grad φ ′ ) Ωc + ∂ ∂t (εgrad φ), grad φ ′ Ω = 0, ∀φ ′ ∈ H 1 Γ E (Ω). (1.86) 1.3. FINITE ELEMENT METHOD
We consider that the electric permittivity does not depend on time and that we can swap the operation grad and ∂ ∂t .

MQS problem in A -φ formulation

In this section, the formulation on A -φ of the magnetic transient (1.51) is studied.

Integrals of the parts by Green's formulas are applied to this formulation to deduce the variational A -φ formulation.

For the first equation of (1.51a), we recall that A magnetic vector potential is defined in H(curl, Ω), and φ scalar potential defined in H 1 (Ω), and let

A ′ ∈ H Γ B (curl, Ω)
be a test function. By applying Green's formula (A.3):

1 µ curl A, curl A ′ Ω + σ ∂A ∂t , A ′ Ωc + (σgrad φ, A ′ ) Ωc = 0. (1.87)
Now, applying the relation of Green (A.2), to the potentials A and φ and to a test

function φ ′ defined in H 1 Γ E (Ω c ), we obtain σ ∂A ∂t , grad φ ′ Ωc + (σgrad φ, grad φ ′ ) Ωc = 0. (1.88)
Finally, the variational form A -φ of the MQS is written as

Find (A, φ) ∈ (H Γ B (curl, Ω) × H 1 Γ E (Ω c )) such as                    1 µ curl A, curl A ′ Ω + σ ∂A ∂t , A ′ Ωc + (σgrad φ, A ′ ) Ωc = 0, σ ∂A ∂t , grad φ ′ Ωc + (σgrad φ, grad φ ′ ) Ωc = 0, ∀(A ′ , φ ′ ) ∈ (H Γ B (curl, Ω) × H 1 Γ E (Ω c )).
(1.89)

Darwin problem in A -φ formulation

Similar to MQS, but it should be mentioned here that the scalar potential φ is defined in the whole domain Ω not just in the conductive domain Ω c , and that only one term is added to each equation of (1.89).

On the same way, for the equation (1.61a), let A ′ ∈ H Γ B (curl , Ω) be a test function.

By applying Green's formula (A.3):

1 µ curl A, curl A ′ Ω + σ ∂A ∂t , A ′ Ωc +(σgrad φ, A ′ ) Ωc + ∂ ∂t (εgrad φ), A ′ Ω = 0.
(1.90) (1.91)

Finally, the variational form A -φ of the Darwin model is written as

Find (A, φ) ∈ (H Γ B (curl, Ω) × H 1 Γ E (Ω)) such as                              1 µ curl A, curl A ′ Ω + σ ∂A ∂t , A ′ Ωc + (σgrad φ, A ′ ) Ωc + ∂ ∂t (εgrad φ), A ′ Ω = 0, σ ∂A ∂t , grad φ ′ Ωc + (σgrad φ, grad φ ′ ) Ωc + ∂ ∂t (εgrad φ), grad φ ′ Ω = 0, ∀(A ′ , φ ′ ) ∈ (H Γ B (curl, Ω) × H 1 Γ E (Ω)).
(1.92)

In the next section, it will be studied how it is possible to approximate such a problem thanks to the finite element method.

Discrete modeling

In the field of computational electromagnetism, Maxwell's equations cannot always be solved analytically, especially in the case when complex geometries are presented.

The numerical approximations of the solutions are widely used via different numerical methods. In the literature, there are several numerical methods to solve the partial differential equations. We cite the finite difference method [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[END_REF], the finite integrate technique (FIT) [START_REF] Weiland | A discretization model for the solution of Maxwell's equations for six-component fields[END_REF], the finite volume method [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF], and the finite element method (FEM) [START_REF] Monk | Finite Element Methods for Maxwell's Equations[END_REF][START_REF] Bossavit | Électromagnétisme, en vue de la modélisation[END_REF], etc.

In the rest of this chapter, we will present the FEM, which allows to pass from an infinite domain to another finite one and to obtain approximated models starting from a spatial and time discretization, in order to apply it to various examples from electrical engineering.

Finite element method

The FEM is one of the most widely used methods in the field of electrical engineering.

Then, we briefly present the concept of this method and the difference between

1.3. FINITE ELEMENT METHOD
the solution of a strong formulation called strong solution and the solution of a weak formulation called weak solution. We consider that such a system of partial differential equations defined in Ω.

Find the solution u of the system defined by

A(u) = f (1.93) Strong Formulation
where A is a differential operator and f is a given continuous function. This formulation mentioned in (1.93) is called the strong formulation.

In order to find the approximated solution of the system (1.93), the use of the numerical methods, in particular the FEM, is necessary. We suppose that we have a Banach functional space V ⊂ Ω in which we look for the exact solution u. In addition, let us take u h in V h , where V h is a discrete functional space included in V .

Applying the variational formulation to (1.93), we have

Find u ∈ V such as a(u, v) = l(v), ∀v ∈ V. (1.94) Weak Formulation
The FEM consists in giving an approximate solution of (1.94). Finally, the discrete problem consists to

Find u h ∈ V h such as a(u h , v h ) = l(v h ), ∀v h ∈ V h (1.95) Discrete Formulation
where a : V × V → R is a bilinear form and l : V → R is linear function. We can assure the existence of a unique solution of the system (1.93) by the Lax-Milgram theorem [START_REF] Monk | Finite Element Methods for Maxwell's Equations[END_REF]. It suffices to show that a is continuous and coercive on V , and L is continuous.

The system (1.95) can be laid as a linear system. The solution

u h is indeed decom- posed into a basis (w 1 , w 2 , ..., w n ) of V h such that u h = n i=1 u i w i (1.96)
with u i the unknowns to be determined and n the degree of freedom (DOF).

As summary, the principle of the finite element method is based on two important steps:

• The first step is used to rewrite the partial differential equations (1.93) in a weak formulation (1.94) showing a bilinear form.

• The second step is the transposition of the system obtained by the first step in a discrete space V h approaching V in order to obtain a system (1.95) allows us to write the problem in matrix formulation.

In the next sub-sections, the functional spaces and discrete variational formulations with the FEM will be studied.

Discrete functional spaces

The use of the FEM leads to a spatial discretization of the domain into geometric elements which can for example be a tetrahedron, a prism or a hexahedron. This spatial discretization is known by the mesh shown in Fig. 1.4. To make computations on meshes, it is necessary to build spaces of discrete functions based on the infinite spaces defined in the part 1.3.1.1. The discrete functional spaces presented in this manuscript were introduced by Hassler Whitney in 1957 [START_REF] Whitney | Geometric Integration Theory[END_REF], in order to be adapted to the discretization of formulations developed previously. This spaces are called the Whitney Elements. We suppose that K is an element of geometry. Let Ω h be the discretization of the domain Ω, produced by the assembly of geometric elements K. The first order discrete functional sub-spaces W 0 (Ω), W 1 (Ω), W 2 (Ω), and W 3 (Ω) are defined as

1.3. FINITE ELEMENT METHOD sub-spaces of H 1 (Ω), H(curl, Ω), H(div, Ω), and L 2 (Ω) as follows W 0 (Ω) = {v h ∈ H 1 (Ω); ∃(n 1 , • • • , n Nn ) ∈ R, v h = Nn k=1 n k w 0 k }, (1.97) W 1 (Ω) = {v h ∈ H(curl, Ω); ∃(e 1 , • • • , e Ne ) ∈ R, v h = Ne k=1 e k w 1 k }, (1.98) W 2 (Ω) = {v h ∈ H(grad, Ω); ∃(f 1 , • • • , f N f ) ∈ R, v h = N f k=1 f k w 2 k }, (1.99) W 3 (Ω) = {v h ∈ L 2 (Ω); ∃(v 1 , • • • , v Nv ) ∈ R, v h = Nv k=1 v k w 3 k } (1.100)
where N n , N e , N f , and N v are the cardinals of the set of nodes, edges, facets, and volumes, respectively. Analogously to the Rham complex defined for continuous functional spaces, there is a complex here called the Whitney complex such as

W 0 (Ω) grad ===⇒ W 1 (Ω) curl = == ⇒ W 2 (Ω) div = = ⇒ W 3 (Ω).
(1.101)

Discrete variational formulation in space

In this part and as we have already introduced the discrete functional spaces, the potential formulations of the quasistatic models given in the section 1.2.2 are studied with the finite element method. We calculate the discrete formulations of weak formulations defined in (1.86), (1.89), and (1.92) in the section 1.3.1.2 in order to make them in matrix form allows to obtain linear systems. All the magnetic and electric quantities are defined in the discrete domains, in the spaces defined in section 1.3.2.2, i.e. they are defined in Ω h . We focus in this section on the Darwin discrete formulation while the discrete forms of the EQS and the MQS models are presented in the appendix C.

Spatial discretization of the Darwin model

Let us introduce φ h the approximated electric scalar potential of φ defined in the space W 0 (Ω h ) ⊂ H 1 (Ω) and A h the approximate magnetic vector potential of A

defined in the space W 1 (Ω h ) ⊂ H(curl, Ω). The weak formulation (1.92) in W 1 (Ω h ) is written Find (A h , φ h ) ∈ W 1 (Ω h ) × W 0 (Ω h ) such as ( 1 µ curl A h , curl A ′ h ) Ω h + (σ ∂A h ∂t , A ′ h ) Ωc h +(σgrad φ h , A ′ h ) Ωc h + ( ∂ ∂t (εgrad φ h ), A ′ h ) Ω h = 0, (σ ∂A h ∂t , grad φ ′ h ) Ωc h + σgrad φ h , grad φ ′ h Ωc h + ( ∂ ∂t (εgrad φ h ), grad φ ′ h ) Ω h = 0. (1.102)
We introduce the decomposition of A and φ in their canonical basis:

A h = Ne i=1 A i w 1 i , (1.103) φ h = Nn i=1 φ i w 0 i . (1.104) Now, we replace φ ′ h by w 0 j a vector of bases in W 0 (Ω h ) and A ′ h by w 1 j a vector of bases in W 1 (Ω h ). The equation (1.102) reads Find (A 1 , ..., A Ne ) × (φ 1 , ..., φ N n ) ∈ R × R such as (1.105)                                      Ne i=1 A i 1 µ curl w 1 i , curl w 1 j Ω h + Ne i=1 ∂A i ∂t σw 1 i , w 1 j Ωc h + Nn i=1 φ i σgrad w 0 i , w 1 j Ωc h + Nn i=1 ∂φ i ∂t εgrad w 0 i , w 1 j Ω h = 0, ∀j = 1, ..., N e Ne i=1 ∂A i ∂t σw 1 i , grad w 0 j Ωc h + Nn i=1 φ i σgrad w 0 i , grad w 0 j Ωc h + Nn i=1 ∂φ i ∂t εgrad w 0 i , grad w 0 j Ω h = 0, ∀j = 1, ..., N n .
The unknowns of the problem are now carried by φ i , i = (1, ..., N n ), and

A i , i = (1, ..., N e ). Let X φ ∈ R Nn and X A ∈ R Ne be the component vectors of φ i and A i ,
respectively. The problem comes down to solving the matrix system

  M C 1 0 G     X A X φ   +   S K C 2 L     ẊA Ẋφ   =   0 0   (1.106)
where ẊA and Ẋφ are the time derivatives of X A and X φ , respectively. In addition, M and S are matrices of size R Ne×Ne ; G and L are matrices of size

R Nn×Nn ; C 1 , C 2 , 1.3. FINITE ELEMENT METHOD
and K are matrices of size R Nn×Ne , R Ne×Nn , and R Nn×Ne respectively defined by:

M i,j = Ω 1 µ curl w 1 i • curl w 1 j dΩ, (1.107) S i,j = Ωc σw 1 i • w 1 j dΩ c , (1.108) G i,j = Ωc σgrad w 0 i • grad w 0 j dΩ c , (1.109) L i,j = Ω εgrad w 0 i • grad w 0 j dΩ, (1.110) C 1i,j = Ωc σgrad w 0 i • w 1 j dΩ c , (1.111) C 2i,j = Ωc σw 1 i • grad w 0 j dΩ c , (1.112) K i,j = Ω εgrad w 0 i • w 1 j dΩ. (1.113)
It should be mentioned here that C 2 is the transpose of C 1 . So, for the sake of simplicity, we will denote them by C and C t respectively.

In the following, the matrix system of the Darwin model is presented in both frequency and time domains.

Discrete variational formulation in frequency domain

In the frequency domain or in harmonic regime, all the quantities depending on time are represented as being the real part of a complex quantity. In addition, the operator of the time derivative ∂ ∂t becomes the factor jω where ω is the angular frequency. Then, the A -φ formulation of the Darwin model is written in frequency domain as

curl 1 µ curl A + σ(jωA + grad φ) + jωεgrad φ = 0, (1.114a) div (σ(jωA + grad φ) + jωεgrad φ) = 0. (1.114b)
Now, applying the finite element method to (1.114) system which implies an asymmetric complex system reads as follows

      M + jωS C t + jωK C 1 jω G + L            X A X φ      =      0 0     
.

(1.115)

Discrete variational formulation in time domain

The system (1.106) is discretized in space but continuous in time. Solving this equation requires the use of numerical temporal integration schemes. We consider

0 = t 0 ≤ t 1 ≤ t 2 ≤ • • • ≤ t i ≤ • • • ≤ t n-1 ≤ t n = T. (1.116)
The time mesh is assumed to be constant:

∀i ∈ {0, ..., n}, ∆t = t i+1 -t i = cst. The time mesh corresponds to a step or time interval ∆t which differs from the stake of time or instant t i defined by t i = i∆t. In this thesis, the time derivatives are approximated by the method of Euler. Euler's method is a numerical method for solving ordinary differential equations. We assume that t ∈ [0, T ], and we have a function y : [0, T ] → R n . Euler's method consists in solving a problem having y a solution of the first order differential equation

     ∂y ∂t = f (t, y(t)), ∀t ∈ [0, T ], y(0) = y 0 (1.117)
where f and y 0 are given. This problem is called Cauchy's problem. Then, by the Cauchy-Lipschitz theorem, we can ensure the existence of a unique solution y to problem (1.117). It suffices to show that the function f is continuous Lipshcitizian uniformly at t. Besides Euler's methods, there are quite a few numerical solving methods. We cite the multiple step methods (midpoint diagram, Adams-Bashforth methods, Adams-Moulton methods and predictor-corrector diagram) [START_REF] Quarteroni | Méthodes Numériques: Algorithm, analyse et applications[END_REF] and single-step methods (the explicit and implicit Euler methods, Runge-Kutta order two and four, etc.). Indeed, the multiple step methods are more precise, but they require more computation steps which presents a drawback

1.3. FINITE ELEMENT METHOD
for the solving of large problem. For that, we are interested in the implicit Euler method. Then, its algorithm is given by

     y i+1 = y i + ∆tf (y i+1 , t i+1 ), y(0) = y 0 (1.118)
with y i = y(t i ).

In the following, only the use of Euler method applied to the Darwin model is presented, those of EQS and MQS are given in the appendix D.

Time discretization of the Darwin model

At time step t i+1 , the derivative of the electric scalar potential φ is written

∂φ ∂t i+1 = φ i+1 -φ i ∆t (1.119)
and the derivative of the magnetic vector potential A is written

∂A ∂t i+1 = A i+1 -A i ∆t .
(1.120)

Applying the implicit Euler scheme (1.118) to (1.106), the system of Darwin model is obtained

       M + 1 ∆t S C t + 1 ∆t K 1 ∆t C G + 1 ∆t L             X A i+1 X φ i+1      = 1 ∆t      S K C L           X A i X φ i     
.

(1.121)

It is interesting to mention here that the resulting Darwin matrix is asymmetric as well as it is ill-conditioned.

The full set Maxwell model [START_REF] Hiptmair | A robust Maxwell formulation for all frequencies[END_REF][START_REF] Badics | Full wave potential formulation with low-frequency stability including ohmic losses[END_REF]] can be invoked to describe the electric and magnetic fields. However, the full set of Maxwell model sometimes introduces the instability of the numerical issue at low frequency [START_REF] Hiptmair | A robust Maxwell formulation for all frequencies[END_REF][START_REF] Ostrowski | H-matrixbased operator preconditioning for full Maxwell at low frequencies[END_REF]. For the A -φ potential formulation proposed in [START_REF] Hiptmair | A robust Maxwell formulation for all frequencies[END_REF][START_REF] Ostrowski | Transient Full Maxwell Computation of Slow Processes[END_REF], the Coulomb gauge (div A = 0) is adopted with Ampere's law and to remove the instabilities in low frequencies. In addition, a generating system approach is used, which is tedious in implementation, and the resulting matrix is singular and difficult to solve. However, the gauged potential formulation proposed in [START_REF] Macneal | A general finite element vector potential formulation of electromagnetics using a time-integrated electric scalar potential[END_REF]] is symmetric but it is only valid for the approximation of the nodal elements of A. Recently, another high frequency gauged potential formulation has been reported in [START_REF] Badics | Full wave potential formulation with low-frequency stability including ohmic losses[END_REF] where edge elements can be used, but the resulting system matrix is asymmetric.

Concerning the formulation A -φ resulting from the Darwin model, due to the very contrasting material coefficients and the magnitude of different types of parameters, the resulting system of singular matrix equations is still ill-conditioned as reported in [START_REF] Koch | A low-frequency approximation to the Maxwell equations simultaneously considering inductive and capacitive phenomena[END_REF]. Moreover, due to the removal of the secondary order terms of A from the original full set of Maxwell equations (1.40), the resulting system matrix is no longer symmetric. In [START_REF] Zhao | A novel gauged potential formulation for 3-D electromagnetic field analysis including both inductive and capacitive effects[END_REF], a novel potential formulation for low frequency applications taking into account the inductive and capacitive effects by adding Coulomb-type gauge as a third equation is presented, which ensures the symmetry of the matrix without additional regularizations. The study is made in both temporal and frequency domains. Recently, a two-step method was proposed in [START_REF] Clemens | A Darwin time domain scheme for the simulation of transient quasistatic electromagnetic fields including resistive, capacitive and inductive effects[END_REF] to avoid the solving of the FE asymmetric and ill-conditioned matrix. The EQS and MQS models are solved iteratively at each time step.

Meanwhile, there is no unique solution of A for the formulation presented in (1.121).

Besides, to ensure the uniqueness of A it would be necessary to add gauge conditions.

Uniqueness of the solution

The uniqueness of the solution of the electric scalar potential φ can be guaranteed either by fixing the average value of φ on the domain Ω such that Ω φ dΩ = 0 or with the boundary conditions imposed on Ω. On the other hand, the uniqueness of the magnetic potential A is not ensured. In order to guarantee the uniqueness of A, from a theoretical point of view, one can impose a Coulomb gauge or from a numerical point of view, a tree gauge [Albanese and Rubinacci, 1990]. These two gauging techniques will be presented in the next sub-sections.

Tree Gauge

The tree gauge or the edge gauge is a gauge acting on the edge elements. The notion of tree comes from the theory of graph. We say in graph theory that a tree is a graph where no path of the graph has the same endpoints and for any pair of vertices in the graph, there exists a path with these vertices for ends (a graph is a set of vertices connected to each other by a set of edges and a path between two vertices of a graph is a serie of edges of the graph, with the two vertices as endpoints of the path). This gauge is defined as follows

A.w = 0 (1.122)
where A is the vector potential and w is an arbitrary vector. Therefore, the tree gauge consists in canceling the component of the magnetic vector potential according to the direction of w. So there is only one path between any two vertices of a tree.

To build a tree from a graph, a graph corresponds to the mesh Ω h , it is necessary that each edge of the graph be associated with a weight. We say here that the mesh or the graph is weighted. Theoretically, the choice of a tree and its weights is arbitrary.

Nevertheless, in [START_REF] Ren | Solving 3D Static Field Problems by Dual Formulations Using Potential Variables[END_REF], it is proved that by solving the matrix problem with a Krylov method, whatever the local minimum towards which this method converges, the magnetic flux density reconstructed from the solution will be unique. In addition, the formulation will be auto-gauged using the conjugate gradient [START_REF] Ren | Influence of the RHS on the convergence behaviour of the curl-curl equation[END_REF].

Coulomb Gauge

The Coulomb gauge is the first possible gauge condition. It means that the normal component of the vector potential A is continuous:

div A = 0. (1.123)
Since the divergence of the edge element basis function is zero within each element, this gauge condition is used when the vector potential is discretized with nodal elements and it is not suitable for edges elements. The large contrast in the values of the material coefficients are encountered and can produce unacceptable errors at material interfaces or geometrical singularities [START_REF] Preis | Different finite element formulations of 3-D magnetostatic fields[END_REF]. Recently, a magnetic vector potential formulation with Coulomb gauge for Full-wave Maxwell problems using edge elements is proposed in [START_REF] Hiptmair | A robust Maxwell formulation for all frequencies[END_REF] and show a large stability and accuracy for all frequencies.

In this thesis, all gauges are used as well as results obtained by ungauged formulations are presented later.

Associated excitation source

It is essential to be able to compute or impose different quantities such as an electric current in a conductor, a difference of electric potentials, a magnetic flux, etc.

In this thesis, we impose an electric potential difference or an electric current in a conductor. In this section, the tools allowing to impose global quantities are developed. To impose a potential difference of a scalar function, one vector β and one scalar α are presented [START_REF] Dular | Coupling of local and global quantities in various finite element formulations and its application to electrostatics, magnetostatics and magnetodynamics[END_REF]]. (1.127)

By the fact of (1.125), we obtain a potential difference of α is equal to 1 between Γ c 1 and Γ c 2 . In fact

γ β.dγ = γ -grad α.dγ = α Γc 1 -α Γc 2 .
(1.128)

The voltage V s is imposed between these two ports, as follows

γ E • dγ = V s (1.129)
and the intensity of the current which enters Γ c 1 [Henneron et al., 2006;[START_REF] Tang | Residual a Posteriori estimator for magnetoharmonic potential formulations with global quantities for the source terms[END_REF] is imposed as follows

Ω J + ∂D ∂t • β dΩ = - Ω J + ∂D ∂t • grad α dΩ = I.
(1.130)

Darwin model

In order to add the terms introduced for the source in the weak formulations, it is necessary to rewrite the electric field E. We decompose E into two terms, E u is the unknown electric field and E s is the source electric field. The latter allows us to impose a voltage in the expression of E.

Voltage

Since γ E s .dl = V s , then β and E s have the same properties to a V s close. Hence, we can write

E s = βV s .
(1.131)

The electric field is written now as follows

E = E u + E s = - ∂A ∂t -grad φ -V s grad α. (1.132)
Then, the strong formulation of the Darwin model with imposed voltage V s reads

           curl 1 µ curl A + σ ∂A ∂t + grad φ + ∂ ∂t (εgrad φ) = -σV s grad α - ∂V s ∂t (εgrad α), div σ ∂A ∂t + grad φ + ∂ ∂t (εgrad φ) = -div σV s grad α + ∂V s ∂t (εgrad α) .
(1.133)

We introduce the decomposition of α in its canonical basis:

α h = Nn i=1 α i w 0 i . (1.134)
Let us denoting X α ∈ R Nn is the component vectors of α i . Consequently, we will have the matrix system of Darwin model (1.121) is written as follows

       M + 1 ∆t S C t + 1 ∆t K 1 ∆t C G + 1 ∆t L             X A i+1 X φ i+1      = -        V s i+1 CX α + ∆V s i+1 ∆t KX α V s i+1 GX α + ∆V s i+1 ∆t LX α        + 1 ∆t      S K C L           X A i X φ i      (1.135)
where ∆V s = V s (t i+1 )-V s (t i ) and M, S, C, K, G, and L are defined in (1.107), (1.108),

(1.111), (1.113), (1.109), and (1.113), respectively, on the page 43.

Current

If we impose an current I in the conductor, the voltage V s would become unknown.

From the equation (1.130), we obtain

Ω σ ∂A ∂t + grad φ + ∂ ∂t (εgrad φ) + σV s grad α + ∂V s ∂t (εgrad α) •grad α dΩ = -I.
(1.136) Finally, the system of equations reads

                     curl 1 µ curl A + σ ∂A ∂t + grad φ + ∂ ∂t (εgrad φ) + σV s grad α + ∂V s ∂t (εgrad α) = 0, div σ ∂A ∂t + grad φ + ∂ ∂t (εgrad φ) + σV s grad α + ∂V s ∂t (εgrad α) = 0, Ω σ ∂A ∂t + grad φ + V s grad α + ∂ ∂t ε(grad φ + V s grad α) • grad α dΩ = -I.
(1.137)

Then, the problem comes down to solve the matrix system

           M C t X t α C t 0 G X t α G t 0 GX α X t α GX α                       X A X φ V s            +            S K X t α K t C L X t α L t CX α LX α X t α LX α                       ẊA Ẋφ Vs            =            0 0 -I           
.

(1.138)

The unknowns of the problem are now carried by (A, φ, V s ).

              M + 1 ∆t S C t + 1 ∆t K X t α (C t + 1 ∆t K t ) 1 ∆t C G + 1 ∆t L X t α (G t + 1 ∆t L t ) 1 ∆t CX α (G + 1 ∆t L)X α X t α (G + 1 ∆t L)X α                          X A i+1 X φ i+1 V s i+1            =            0 0 -I i+1            + 1 ∆t            S K X t α K t C L X t α L t CX α LX α X t α LX α                       X A i X φ i V s i           
.

(1.139)

1.5. NUMERICAL VALIDATION
The two matrix systems with different excitation sources using the Darwin model have been presented in this part while it is presented for the EQS and MQS formulations in the appendix E.

Numerical validation

Model presentation

Geometry

An application concerning a conductive cylindrical geometric shape is considered in Fig. 1.7 which is made of aluminum, having a section of 3.14 mm 2 and a length of 100 mm, and an air box encloses the electric cable. To model the capacitive effects in this application, the electric potential on the below face of the air box, marked on blue, is taken to be zero which corresponds to connect to the ground.

Material characteristics

The studied geometry is composed of two parts, a conductive part representing the electric cable in red part and the non-conductive air box in grey part as shown in Fig. 1.7. It should be mentioned here that the electric permittivity of the air has been increased in order to approach the appearance of capacitive effects at reasonable frequencies. All the characteristics for these materials are presented below:

1) The electromagnetic device is encompassed by an air box having a:

• relative magnetic permeability µ r = 1, Figure 1.7: The geometry of the electric cable.

• relative electrical permittivity ε r = 1000.

2) The aluminum having a:

• relative magnetic permeability µ r = 1,

• relative electric permittivity ε r = 1,

• electric conductivity σ = 20 MS/m.

Computational configurations

To take into account the skin effects in the conductor for high frequencies, the geom-

Evolution of the impedance versus the frequency

The evolution of the impedance as a function of the frequency obtained from EQS, MQS, and Darwin models as well as the phase are presented in Fig. 1.12. For all 1.5. NUMERICAL VALIDATION models, the impedance Z corresponds to the DC resistance of the conductor with Z = R DC = 0.32 Ω when the frequency tends to zero.

Figure 1.12: Modulus of impedance and the phase computed with respect to frequencies.

When the frequency increases, due to the skin effect appearing in the conductor, the value of the resistance of the conductor increases. This effect is not taken into account with the EQS model. For f ∈ [0 : 10 7 ] Hz, the MQS and Darwin models give a similar evolution of the impedance. Indeed, the capacitive effect is negligible compared with the inductive and resistive effects. For f > 10 MHz, the capacitive effect appears and it is responsible for the drop on the evolution of the impedance for both EQS and Darwin models. As expected, the resistive-capacitive coupled effects are presented by the EQS model, while the coupled resistive-inductive effects are modeled by the MQS model. As long as by invoking the Darwin model, the resistive-inductive-capacitive effects can be captured for f > 10 MHz. Then, the resonant frequency of this example is around f = 100 MHz.

Post-processing indicator

Based on the considered example in Fig. 1.7, carried out in order to identify the various quasistatic domains, we compute in post-processing the following electrical powers by using the Darwin model:

P 1 = 1 2 ε||jωgrad φ|| L 2 (Ω)
(1.140)

P 2 = 1 2 ε||ω 2 A + jωgrad φ|| L 2 (Ω) (1.141)
where the P 2 contains the neglected term of the full set Maxwell system which is responsible of the radiation effects. The aim of this part is to specify the frequency for which the neglected term should be taken into account and to investigate the limit of the Darwin model in the frequency domain.

The electrical powers mentioned-above are presented in Fig. 1.13a, which shows that from f = 1 GHz, it is essential to take into account the term modeling the wave propagation effects. In addition, the relative error between the two terms becomes important at this frequency as shown in Fig. 1.13b.

Consequently, based on the parameters defined in (1.36) and the results obtained by the previous application, a diagram specifies the various electromagnetic domains namely, static, quasistatic, and Full-wave as shown in Fig. 1.14. Indeed, when k 1 >> 1, the electric and magnetic fields are decoupled which means that we are in the static regime. In this zone, the electrokinetic, the electrostatic, and the magnetostatic models can be used to compute the resistance R, the capacitance C, and the inductance L, respectively.

On the other hand, the k 2 and k 3 indicators which determine the necessity to use the MQS and EQS models, show an almost parallel variation, indicating that when k 2 moves away from 1, the MQS model cannot be adopted, and when k 3 approaches to 1, the model EQS can be adopted. In particular, these indicators show that when the frequency increases, the solution obtained by MQS prones to errors since the capacitive effects should be modeled.

Besides, the Darwin model can be invoked instead of MQS model when k 1 ≳ 1. In addition, around 3 GHz, the Darwin model reaches its limit. In this case, the study

Conclusion

In this chapter, the Maxwell equations and the material laws as well as the boundary conditions are recalled. Secondly, the different quasistatic models, namely EQS, MQS, and Darwin models are introduced, as well as their potential formulations.

Thirdly, the finite element method associated is presented. Finally, in the last part,

Solvers

Introduction

The solving of linear systems plays an important role in the numerical simulation of many scientific and industrial problems. In this domain, the use of finite element computation codes often leads to sparse linear systems that can reach sizes of tens of millions of unknowns, as the matrix system defined in section 1.3 for EQS, MQS, and Darwin models. The solving of these systems has the disadvantage of leading to prohibitive costs in CPU time and memory space.

Partial differential equations (PDE) appear in several scientific fields: fluid mechanics by Navier Stokes equations, electromagnetism by Maxwell's equations, etc. Indeed, in this domain, numerical simulation plays a primordial role. It aims at determining the electromagnetic fields in response to an excitation source (charge density, current density, etc.). This task is carried out by various techniques and modeling methods which make it possible to solve electromagnetic problems, in the time or frequency domain.

These systems can be symmetric or asymmetric, real or complex, with a single or multiple second member.

In the previous chapter, the A -φ formulations are presented in matrix form, using the FEM for the EQS, MQS, and Darwin models. In general, at each time step, we solve a system of the form

AX = F (2.1)
where A represents the FE matrix defined in the section 1.3.2, X is the vector of unknowns which represents to X φ for the EQS model and (X A , X φ ) for the MQS and Darwin models.

By relying on the FEM, the unknowns X A are the circulations of the magnetic vector potential A along each edge e k , k = 1, • • • , N e , of the mesh Ω h and the X φ are the electric scalar potential φ defined on each node

n k , k = 1, • • • , N n , of Ω h .
To determine these potentials, the system (2.1) should be solved. There are various methods of solving, one can refer to the direct methods or the iterative methods [START_REF] Kelley | Iterative Methods for Linear and Nonlinear Equations[END_REF].

Direct methods

The principle of the direct methods to determine the solution X is not to seek the inversion of the matrix A in order to obtain the solution X = A -1 F because the 2.1. SOLVERS inversion of the matrix A is expensive in time and memory. Then, to determine the solution X, direct solvers exist such as the Cholesky method, the LU decomposition, or the Gauss-Jordan method. The Cholesky method is applicable only for a positive definite symmetric square matrix A. This method computes a lower triangular matrix L such that A = LL t . The LU decomposition is based on a factorization of the matrix A defined in (2.1) into two matrices upper triangular U and lower triangular L. This principle is called Gaussian elimination. To find the solution X, it is then sufficient to first solve the system LY = F and secondly the system UX = Y. The Gauss-Jordan method consists in performing identical operations on the matrix A and the vector F so that the matrix A becomes the identity matrix.

On the other hand, these direct methods can fail for various reasons. For example, if the matrix is very large. Here we are talking about the case with 40-50 million unknowns. Moreover, when the matrix contains a lot of zeros such as the sparse matrices. The main concern of direct methods in this case occurs during the factorization process: operations between non-zero elements generate new ones in the factorized matrix, where there was a zero before, which means that certain terms in the matrices L and U are non-zero while the initial terms of the matrix A are zero.

The factors are more filled than the initial matrix, which is called the phenomenon of filling.

It is now possible to solve linear equations of very large sizes thanks to the development of parallel sparse direct solvers like the MUMPS (MUltifrontal Massively Parallel sparse direct Solver) [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF], PARDISO (PARallel DIrect SOlver) [START_REF] Schenk | Solving Unsymmetric Sparse Systems of Linear Equations with PARDISO[END_REF], and SuperLU [START_REF] Li | SuperLU_DIST: A scalable distributedmemory sparse direct solver for unsymmetric linear systems[END_REF] solvers.

The solutions obtained by these solvers are generally very accurate. We consider in the following the solver Mumps as direct solver to solve the different quasistatic formulations defined in the chapter I.

Iterative methods

Direct methods are preferred since they are more efficient than iterative methods but when the solving of such system is out of support of the direct methods, the iterative solvers are used. The principle of this solvers is to start with an initial solution and then iteratively move towards the correct solution for X. The iterative solver continues to iterate until it finds a solution that satisfies the stop criterion.

It is not an exact solution but iterative solvers are suitable for large computational problems and have rapidly gained popularity in the various fields of scientific computation thanks to their low computational and memory costs, their simplicity of implementation and their adaptability to high performance computation.

Two main classes of iterative solvers exist to solve the system (2.1), one names the iterative solvers of the fixed type and the Krylov sub-space solvers. The fixed point type solvers are simple to implement but they are only provided for certain types of matrices. The Jacobi method and the Gauss-Seidel method are iterative fixed point solvers. Krylov space-based solvers are the most well-known and used solvers. They are widely used for solving sparse linear systems of large sizes [START_REF] Saad | Iterative solution of linear systems in the 20th century[END_REF]. These methods are based on a projection technique on a Krylov subspace, of dimension smaller than the size of the problem. They approach the solution by a vector in a Krylov sub-space. The Conjugate Gradient Solver (CG) is the best known Krylov iterative solver. In addition, there are other iterative Krylov solvers such as BiCG (BiConjugate Gradient method), BiCGSTAB (BiConjugate Gradient STABilized method), and GMRES (Generalized Minimal RESidual), etc.

When the studied matrix is symmetric and positive definite, the CG method [START_REF] Hestenes | Methods of Conjugate Gradients for Solving Linear Systems[END_REF] is the most prominent iterative method for solving sparse systems When the studied matrix is asymmetric, BiCG method is the most used among different Krylov methods. This method was developed by Felcher in 1974. It is a generalization of the CG method. The latter method is the fastest iterative solver, but only for positively defined symmetric systems. Suppose that the Krylov subspace

K n = span {F, AF, • • • , A n F} . (2.
2)

The CG methods look for approximate the solutions at each step n in the Krylov sub-space. In general, the iterative methods allow to achieve the solution of a linear system as the limit of a series of approximate solutions. Unfortunately, this may fail if you apply it excessively. Since the GMRES algorithm minimizes the residue over Krylov spaces. Therefore, by performing a step of the GMRES algorithm after each BiCG step, the resulting iteration is stable; this is generally called BiCGSTAB [van der Vorst, 1992]. Thus, the latest is in principle a more general solver than CG but suffers from lower efficiency when applied to the problems for which CG was intended.

There is one important thing to note here that, the convergence of the iterative solvers is sensitive and there exist a matrix which will cause them to fail to converge.

DIFFERENT FORMULATIONS/CASES

Besides, the finding of a good preconditioner is necessary to accelerate the speed of convergence.

The resulting finite element matrix from the Darwin model is ill-conditioned. Therefore, the speed of convergence of iterative solvers slows down. On the other hand, to improve the conditioning of the matrix, and consequently to accelerate the speed of convergence, it is essential to add effective preconditioners. Among the methods which are used to precondition the linear system, we can cite the relaxation methods, the incomplete factorization methods or the multigrid methods [START_REF] Falgout | The Design and Implementation of hypre, a Library of Parallel High Performance Preconditioners[END_REF].

In this thesis, the preconditioner used is the Split-Jacobi, which corresponds to a relaxation method. The deriving preconditioner splits the Jacobi preconditioner M (M is the diagonal of A) as M = P t P, where P i,j = M i,j . Then, the final system to solve is given as follows

P -1 A(P t ) -1 (P t X) = P -1 F. (2.3)
Therefore, solving (2.1) is replaced by solving (2.3). In the case of Jacobi preconditioner, the system is written M -1 AX = M -1 F. Here we say that the system is left preconditioned. Consequently, the system (2.3) is left and right preconditioned at the same time. The Jacobi and the Split-Jacobi are two completely different preconditioners. The Split-Jacobi is a preconditioner that keeps the symmetry of the matrix while Jacobi breaks it. Consequently, it modifies the matrix more deeply and generally that allows to have a better conditioning of the resulting finite element matrix and, a fortiori, a better convergence.

Different formulations/cases

In the previous chapter, the electromagnetic formulations are studied in both time and frequency domains. In this section, different formulations of the Darwin model are presented in the frequency domain with different cases.

There are three different formulations for the Darwin model will be studied in the following.

M1: first formulation

The classical A -φ formulation of the Darwin model given in [START_REF] Koch | A low-frequency approximation to the Maxwell equations simultaneously considering inductive and capacitive phenomena[END_REF] without entering the gauge condition and with imposed voltage V s is written

       curl 1 µ curl A + σ(jωA + grad φ) + jωεgrad φ = -σV s grad α -jωV s (εgrad α),
div (σ(jωA + grad φ) + jωεgrad φ) = -div (σV s (grad α) + jωV s (εgrad α)) .

(2.4)

Now, applying the finite element method, the MVP A is discretized with edge elements and the ESP φ with nodal elements which introduce an asymmetric matrix.

Then, the matrix system of the ungauged Darwin formulation model in frequency domain reads

     M + jωS C t + jωK jωC G + jωL           X A X φ      = -      V s CX α + jωV s KX α V s GX α + jωV s LX α      =      F A (ω) F φ (ω)      (2.5)
where M, S, C, K, G, and L are the same matrices defined in the chapter I on the page 43.

The system (2.5) is not symmetric, ill-conditionned, and ungauged. Since no gauge condition is applied here, the system does not have a unique solution in terms of the discrete values of the potentials. Consequently, the system is singular. The system is not symmetric due to the dielectric contribution in the upper right block, the counterpart of which is absent in the lower left block. Therefore, iterative solvers suitable for the solution of complex asymmetric systems, such as the BICGSTAB or the GMRES method, should be applied.

In this thesis, the BiCGSTAB solver is used as an iterative solver to solve the nonsymmetric Hermitian system. However, due to the large differences in the values of the material coefficients, the severe bad conditioning could considerably reduce the convergence rate of the iterative solver. Consequently, an adapted preconditioner is needed. Accordingly, among the preconditioners existing in the literature, one applies the method of Split-Jacobi to solve the system in a diagram of preconditioning by blocks.

M2: second formulation

The formulation mentioned above is ungauged. As the right hand side (RHS) of (2.5) is divergence-free per construction, the use of iterative solvers is efficient. In the

DIFFERENT FORMULATIONS/CASES

following, a Coulomb gauged formulation for the Darwin model will be studied. This formulation is proposed in [START_REF] Zhao | A novel gauged potential formulation for 3-D electromagnetic field analysis including both inductive and capacitive effects[END_REF] which is also applied for the fieldcircuit coupled formulation [START_REF] Zhao | A novel gauged potential formulation for 3-D electromagnetic field analysis including both inductive and capacitive effects[END_REF]. By introducing the Lagrange multiplier p, and rewriting the current continuity (second equation) equation using the gauge technique, the resulting matrix becomes symmetric. Then, the Coulombtype gauge A -φ formulation for the Darwin model with imposed voltage V s reads

                           curl 1 µ curl A + σ(jωA + grad φ) + jωεgrad φ -jωεgrad p = -σV s grad α -jωV s (εgrad α), div (σ(jωA + grad φ) + jωεgrad φ)+jωdiv (εA) = -div (σV s (grad α) + jωV s (εgrad α)) , jωdiv (εA) = 0. (2.6)
It should be mentioned here that by applying the finite element method, the Lagrange multiplier p is approximated with the nodal elements. Then, the matrix system reads

            M + jωS C t + jωK t -jωN t C + jωK 1 jω G + L 0 -jωN 0 0                        X A X φ X p            =            F A (ω) F φ (ω) 0           
.

(2.7)

The unknowns of the problem are now carried by (A, φ, p) where A and φ are the potentials magnetic and electric unknowns respectively, and p i , i = (1, ..., N n ), is a scalar nodal unknowns. So, let X p ∈ R Nn , the component vector of p. And, N is the same matrix K of size R Nn×Ne defined in (1.113) on the page 43.

The Coulomb-type gauge is additionally imposed as a third equation by the Lagrange multiplier which ensures the symmetry of the matrix without additional regularizations. This fact has a great effect on the solving using the direct methods, where one can greatly reduce the computational time as well as the memory cost.

In [START_REF] Zhao | A novel gauged potential formulation for 3-D electromagnetic field analysis including both inductive and capacitive effects[END_REF], it is noted that it is now convenient and quick to review solvers in the literature thanks to the modern multi-core processor and the drop in memory prices. For example, the solving of a complex and symmetric system with 1.38 million unknowns by PARDISO [START_REF] Schenk | On fast factorization pivoting methods for sparse symmetric indefinite systems[END_REF]] takes around 201 s using two Intel Xeon processors of type E5-2687W v3 (10 cores, 25 MB cache, 3.10 GHz) and 38 GB of memory (approximately 28 GB of memory for storing the factors of the finite element matrix) additional cost as indicated in [START_REF] Zhao | A new stable full-wave Maxwell solver for all frequencies[END_REF]. On the other hand, the cost of the memory is high at 65 GB of memory (about 56 GB of memory to store the factors) and the computational time at about 440 s using a complex asymmetric PARDISO solver. To conclude, the symmetric gauged formulations are preferable because during the storage of the factors of the finite element matrix, one is able to save approximately 50% of memory as well as the execution of numerical factorizations is much faster.

However, in the case of industrial applications, the number of DoFs is huge and in practice, it is difficult to use the direct solvers due to the memory limitation.

Thus in our work, we choose iterative solvers with an efficient preconditioner. They are less memory intensive than direct methods because they only use vector matrix products. For large systems, when properly chosen and configured, they are also more computationally efficient. Nevertheless, the linear system (2.7) is a classical saddle point problem which leads to a system of the form

  A 1 B t B 0   A   x y   =   F 1 F 2   (2.8)
where A 1 is a matrix of size n × n and B is of dimension m × n and F 1 and F 2 are the RHS vectors.

Although this matrix is not singular, solving the system by iterative methods using classical preconditionings [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] does not always offer satisfactory convergence rates, sometimes going so far as to diverge. Although the system (2.8) is linear, solving it poses a significant challenge. Indeed, the matrix A has both positive and negative eigenvalues.

M3: third formulation

As mentionned before, this kind of problem (2.8) is found in a wide variety of science and engineering applications. There is a situation leading to a non-zero block is the discretization of the equations describing slightly compressible fluids or solids. In addition, we can cite techniques used in the problems of finance [START_REF] Perold | Sparsity and piecewise linearity in large portfolio optimization problems[END_REF], electromagnetism [START_REF] Hallquist | Sliding interfaces with contactimpact in large-scale Lagrangian computations[END_REF], optimization under constraints [START_REF] Wright | Primal-Dual Interior-Point Methods[END_REF], and the Mortar method [START_REF] Rodger | Coupled elements for problems involving movement (switched reluctance motor)[END_REF][START_REF] Bouillault | The Mortar edge element method in three dimensions: Application to magnetostatics[END_REF][START_REF] Wohlmuth | A Mortar finite element method using dual spaces for the Lagrange multiplier[END_REF][START_REF] Antunes | Using hierarchic interpolation with Mortar element method for electrical machines analysis[END_REF]]. What's more, a review of the literature on methods of solving saddle point problems can be found in Benzi and al. (see [START_REF] Benzi | Numerical solution of saddle point problems[END_REF]).

The system to be solved is a saddle point problem, which requires ad-hoc solvers and

DIFFERENT FORMULATIONS/CASES

preconditioners [START_REF] Benzi | Numerical solution of saddle point problems[END_REF]. Whereas, one seeks to adapt the M2 formulation in order to be able to apply the iterative solvers.

The method of Lagrange double multipliers [START_REF] Aubertin | Single and double Lagrange multipliers approaches applied to scalar BIBLIOGRAPHY potential formulation in magnetostatic fem[END_REF] can be applied in our case in order to improve the conditioning of the matrix. This method consists in duplicating y such that: y = y 1 + y 2 and with y 1 = y 2 . Then, the system to be solved reads in the following form

     A 1 B t 0 B I -I 0 -I I      A      x y 1 y 2      =      F 1 F 2 F 3     
.

(2.9)

Thanks to this method, the system (2.8) can be solved now using iterative methods when direct solvers do not meet our needs.

Let us return to the Darwin model especially to the formulation (2.7), and by applying the method of double multipliers of Lagrange, the modified symmetric Coulomb gauge formulation with imposed voltage V s reads

                  M + jωS C t + jωK t -jωN t -jωN t C + jωK 1 jω G + L 0 0 -jωN 0 I -I -jωN 0 -I I                                    X A X φ X p 1 X p 2                  =                  F A (ω) F φ (ω) 0 0                 
.

(2.10)

The unknowns of the problem are carried by (A, φ, p 1 , p 2 ) where p j i , i = (1, ..., N n ), j = 1, 2 are the scalar nodal unknowns. And, X p j ∈ R Nn are the components vector of p j .

Summary

To Moreover, the tree gauge can be adopted on the M1 formulation to ensure the uniqueness of A (see the section 1.3.3). In the following, seven cases are taken into account as shown in The solving with iterative solvers without adding a preconditioner may sucessfully have an approximate solution, but often it is prone to errors when the matrix system is ill-conditioned, in particular, in the case of the Darwin formulation. However, by adding the Split-Jacobi preconditioner, the BiCGSTAB solver converges when applied for the asymmetric and ungauged M1 formulation, as well as in the case of incorporating the tree gauge. On the other hand, the Mumps solver can be adapted for the asymmetric M1 formulation and for the Coulomb gauged M2 formulation where the resulting system is symmetric. As long as the M2 formulation can only be solved by Mumps, the M3 formulation is proposed with an iterative solver in the case of the large size problem. In the following part, the results obtained by each case will be studied and compared.

Comparison of formulations with different solvers in the frequency domain

Numerical results are presented for a parallel plate capacitor as shown in Fig. 2.1a

with alternating voltage excitation at the plate terminals considering all the cases mentioned before. The used mesh, presented in Fig. 2.1b, features 57,237 elements including 9,879 nodes and 67,474 edges. For all the test cases, the magnitude of the voltage excitation is taken to be V smax = 1 V. The frequency interval is [1 : 10 9 ] Hz.

The electric conductivity of the conductor is taken as 50 kS/m. For the relative permittivity, it is set as 1 for the conductor (grey part) and 10 for the dielectric (green part), while the relative permeability is set as 10 3 and 1, respectively. In Table 2.2, the DoFs of the M1, M2, and M3 formulations are presented according to the seven cases proposed above. The DoFs vary depending on the different cases such as the Coulomb gauge and the tree gauge as well as the addition of a Lagrange multiplicator which will duplicate the number of nodal unknowns.

In the following, the current flowing out of the terminal of the capacitor is computed using the M1, M2, and M3 formulations. In [START_REF] Zhao | Robust full-wave Maxwell solver in time-domain using magnetic vector potential with edge elements[END_REF], it is shown that the current stability of all methods in this interval. When f > 10 MHz, a difference between the results can be observed as shown in Fig. 2.2b, but all results remain close to the reference. We notice that the M1 formulation is the least stable among the others.

The ungauged formulation solved with the BiCGSTAB is far from the reference as well as that gauged with the tree gauge and solved by Mumps.

On the other hand, the M1 formulation gauged by the tree gauge presents good results when it is solved by the BiCGSTAB solver, but it is too expensive in computational time as shown in Fig. 2.3. Meanwhile, a good behavior is presented by the M2 formulation when the Mumps solver is adopted. It should be mentioned here that the M1 formulation solved with Mumps presents identical results as the one obtained by M2, but the solving of the M1 formulation is more expensive in memory and computational time since the FE matrix is asymmetric.

Besides, the M3 formulation, the modified symmetric formulation, presents good results by comparing to the ungauged formulation when the solving is done with iterative methods. This formulation shows an improvement in results by comparing with the results obtained by the M1 formulation. In conclusion, among the cases cited in Table 2.1, the gauged and symmetric formulations are ideal whether solving with direct or iterative methods, more precisely, the M2 and M3 formulations.

The computational time with direct solver, namely Mumps, is stable when the frequency increases. On the other hand, the computational time for BiCGSTAB iterative solvers varies, which depends on frequency and the conditioning number of the matrix. The use of non-symmetric M1 formulation with a tree gauge is the most expensive as shown in Fig. 2.3. In general, the direct solver Mumps is faster than the iterative solver BiCGSTAB especially for problems with several millions unknowns and if the system is symmetric, we can save 50% of computational time.

In addition, the M3 formulation needs much more time than M1 formulation because the number of nodal unknowns is three times bigger as shown in Table 2.2. The number of iterations and the residue norm for the iterative solver, namely BiCGSTAB with Split-Jacobi preconditionner, are presented in the following with the different formulations.

Performance of M3 formulation

In this part, we study the performance of the M3 formulation, i.e. the convergence of the BiCGSTAB iterative solver in comparison with the different cases resulting from the M1 formulation as well as the number of iterations necessary to obtain the solution.

The number of iterations required by each case as a function of frequency is shown in Fig. 2.4. It is clear that the non-symmetric formulation gauged by the tree technique requires a large number of iterations to obtain an acceptable solution, this shows that although this case offers us acceptable solutions but it is too expensive since by applying this technique the resulting finite element matrix is more ill-conditioned.

On the other hand, for the two other cases, the ungauged asymmetric M1 formulation and the symmetric M3 formulation with the Coulomb gauge, require a close number of iterations, as well as much faster than the other case.

Furthermore, we choose two frequencies in order to compare the behavior of the But, the two other propositions present a super linear behavior and a solution with a reduced number of iterations while benefiting from a greater robustness. In addition, the stability of the solution using the M3 formulation is retained for a wide frequency range in comparing with the M1 formulation. The M3 formulation seems well suited for industrial cases where iterative methods are used.

Discussion

For coarse meshes, the direct methods are much faster than the iterative methods due to the fact that the factorization happens quite quickly. When the number of DoFs increases, the LU method becomes more and more expensive, resulting in a very high computational time compared to the other approaches. One thus finds the standard behavior of the direct solvers. Furthermore, the iterative solvers are preferred to solve the Darwin formulation when the resulting matrix is ungauged (2.5). The iterative solvers show a good performance for intermediate frequencies. When using a tree type gauge, both iterative and direct solvers have good results. But looking at Split-Jacobi preconditioner.

COMPARISON OF SOLVERS VERSUS NUMBER OF ELEMENTS

the computational time in Fig. 2.3 and the residue norm in Fig. 2.5, the use of this kind of gauge is too expensive. Alternatively, the symmetric formulation presented in (2.7) is ideal when using direct solver, namely Mumps, but due to the limitation of the direct solver in terms of DoFs, the use of the direct solvers is not guaranteed for industrial applications. Therefore, we resort to the iterative solver when the matrix is symmetric by using the M3 formulation proposed in (2.10).

Comparison of solvers versus number of elements

In this part, we study the performance of the solvers as function of the number of elements of the mesh. Although the stability and the performance of the iterative solvers depend on the frequency and the conditioning of the matrix. This study is made out on the capacitor model presented in Fig. 2.1 with the different formulations studied in the previous section where they are solved by the direct and iterative methods. In the following, nine different meshes are considered going from 50 thousand elements up to 15 million elements. This study is done in both frequency and time domains.

Frequency domain

For different frequencies as well as with a variation in number of the elements, the direct solver Mumps will give stable and accurate results. In contrast, the iterative solver depends strongly on the frequency and number of elements. For this reason and according to the results obtained in the previous sections, the results obtained by Mumps will reference for the following study.

To study the performance of the iterative solvers with the M1 and M3 formulations, the modulus of the current flowing out of the terminal of the conductor are presented in Table 2.3. Note that for all frequencies, the results obtained with a direct solver (Mumps) are stable. For f < 10 MHz, the results obtained by BiCGSTAB with a Split-Jacobi preconditioner are accurate and consistent with those obtained by Mumps as shown in Table 2.3a-2.3c. On the other hand, for f = 10 MHz, the results obtained with the iterative solver show an instability when the number of elements increases as shown in Table 2.3d, which is equivalent to a convergence towards the approximate solution was not reached by the BiCGSTAB.

To model the capacitive-inductive effects coupled under the condition of neglecting the wave propagation, the frequencies reached by this study are sufficient, and de- Table 2.3: The performance of iterative solver (BiCGSTAB-M3) with a Split-Jacobi preconditioner and direct solver (Mumps-M2) versus the number of elements for different frequencies.

COMPARISON OF SOLVERS VERSUS NUMBER OF ELEMENTS

pend on each application. So at high frequencies, higher than 10 MHz, the solver mentioned with the Split-Jacobi adopted, can leads to unexpected errors. Therefore, the construction of generally effective and efficient preconditioners is increasingly desirable.

Time domain

In this part, we briefly present the results obtained for the M1, M2, and M3 formulations in the time domain. Since the M1 formulation is already presented in (1.61) in

Chapter I, we begin by introducing the Coulomb-type gauge Darwin formulations.

Then, the M2 formulation in time domain reads

                                         curl 1 µ curl A + σ ∂A ∂t + grad φ + ∂ ∂t (εgrad φ) - ∂ ∂t (εgrad p) = -σV s grad α - ∂V s ∂t (εgrad α), div σ ∂A ∂t + grad φ + ∂ ∂t (εgrad φ) +div ε ∂A ∂t = -div σV s grad α + ∂V s ∂t (εgrad α) , div ε ∂A ∂t = 0.
(2.11)

Then, the matrix system of this formulation in time domain reads

     M C 1 0 0 G 0 0 0 0           X A X φ X p      +      S K N C 2 + K L 0 N 0 0           ẊA Ẋφ Ẋp      = -      V s CX α + Vs KX α V s GX α + Vs LX α 0     
. (2.12)

Applying the implicit Euler scheme (1.118) to (2.12), the matrix system of the Coulomb-type gauge M2 formulation is obtained

              M + 1 ∆t S C t + 1 ∆t K t 1 ∆t N t C + 1 ∆t K ∆tG + L 0 1 ∆t N 0 0                          X A i+1 X φ i+1 X p i+1            = -              V s i+1 CX α + ∆V s i+1 ∆t KX α V s i+1 GX α + ∆V s i+1 ∆t LX α 0              + 1 ∆t            S K t N t C + K L 0 N 0 0                       X A i X φ i X p i           
.

(2.13)

The resulting system matrix (2.13) is symmetric. For the same reasons as in the frequency domain, one introduces the M3 formulation in the time domain which is written solvers will be studied in Table 2.4. In this table, the current is presented for the first time step. Since the direct solver presents always the good results for different frequencies, it is also adapted in the time domain as a reference for the comparison.

                     M + 1 ∆t S C t + 1 ∆t K t 1 ∆t N t 1 ∆t N t C + 1 ∆t K ∆tG + L 0 0 1 ∆t N 0 I -I 1 ∆t N 0 -I I                                       X A i+1 X φ i+1 X p 1i+1 X p 2i+1                  = -                    V s i+1 CX α + ∆V s i+1 ∆t KX α V s i+1 GX α + ∆V s i+1 ∆t LX α 0 0                    + 1 ∆t                  S K t N t N t C + K L 0 0 N 0 I -I N 0 -I I                                   X A i X φ i X p 1i X p 2i                  . ( 2
The stability of the M1 formulation can be obtained by using the tree gauge but the latter technique is too expensive in the numerical term. As shown in Table 2.4, the ungauged M1 formulation exhibits a deviation when comparing to the results obtained by Mumps. Besides, in the case of M3 formulation, where the resulting matrix is symmetric and positive definite, the BiCGSTAB can be reduced to CG, which will bring us an additional speed in the rate of convergence. Then, the stability of the solution is obtained for a frequency of 10 MHz by adapting the M3 formulation instead of M1 one.

It should be mentioned here that in the case where the number of elements exceeds 30 million, the direct solver failed in the factorization step due to the large size of the resulting finite element matrix. Besides, an accurate result is obtained by using the iterative solver with the M3 formulation.

Table 2.4: The performance of iterative solver (BiCGSTAB-M1)-reference and (CG-M3) with a Split-Jacobi preconditioner and direct solver (Mumps-M2) versus the number of elements for f = 10 MHz.

Moreover, in order to check the results for several time steps, two periods of currents

Conclusion

In this chapter, two large families of solvers, direct methods and iterative methods, have been presented. They play the most important role in solving the partial differential equations, in particular, the potential formulations resulting from the electromagnetic problems. In addition, three Darwin formulations are studied, as detailed in section 2.2, M1, M2, and M3 formulations. The numerical results obtained with these different formulations are investigated. For coarse meshes up to a several million, direct solvers are ideal and preferred when the resulting matrix is symmetric, in particular, the case of the M2 formulation. In addition, this formulation is only suitable for the direct methods. On the other hand, in this thesis, we are interested in the modeling of industrial applications at high frequencies, where the number of elements may expand and attract more than 30 million elements. In this situation, the direct solvers are prone to problems since it is out of their support.

Therefore, the M3 Coulomb-type gauge and symmetric formulation is proposed to take advantage of the solving with iterative solvers.

Surge arrester

Introduction

Surge arresters are one of the devices used on the electric energy transmission and the distribution networks to limit the high voltages that occur mainly during lightning strikes. This type of application is arranged in the "dielectric" category. The field effects in such high-voltage devices are capacitive and non-linear resistive effects [START_REF] Hinrichsen | Grundlagen der Feldsteuerung, presented at the Feldsteuernde Isoliersysteme: Werkstoffe, Design, Prüfung und Simulation[END_REF]. The constitutive laws of materials are only described by the electric permittivity and conductivity. Then, they can be adequately described using the EQS model [Haus and Melcher, 1989b].

Surge arresters can be classified into two main families: silicon carbide varistors (SiC) and zinc oxide (ZnO) varistors. In our study, we will concentrate on the ZnO surge arresters since most studies focus on this family [START_REF] Clark | High-power electrical pulses using metal oxide varistors[END_REF][START_REF] Weida | Design of zno microvaristor material stress-cone for cable accessories[END_REF][START_REF] Abd-Rahman | Stress control on polymeric outdoor insulators using zinc oxide microvaristor composites[END_REF]. ZnO surge arresters are used on transmission and distribution networks to protect equipment against overvoltages and to improve service quality. These devices have a strongly non-linear voltage-current V (I) characteristic: they are practically insulating at low voltages and become more and more conductive with increasing voltage at their terminals. Then, to characterize the voltage-current and to capture the capacitive-resistive effects simultaneously of the arresters, the EQS model should be considered. In the literature, the significant 3.2. SURGE ARRESTER evolution over time is reported when a variation of the frequency presents. The capacitive and resistive effects from low to high frequency can be well observed.

The present surge arrester is studied in time and frequency domains. In the time domain, two different cases are considered: linear and non-linear cases. For the linear case, the model is considered in a steady state at 50 Hz. Then, we check that a weak current is flowing through the central ZnO part. Otherwise, we introduce a nonlinear electric constitutive law on the conductivity of ZnO in order to characterize the curve V (I). In the frequency domain, we compute the impedance versus frequency in order to verify the capacitive behavior of the surge arresters.

Model presentation

Geometry

A ZnO surge arrester is made up of a set of cylindrical pellets staked on top of each

Material characteristics

This device can be represented schematically, as shown in Fig. 3.3, by two conductive parts made of aluminum and cast iron which are very good electric conductors, dielectrics, and an air tube. Thereafter, all the characteristics are presented below:

1) The surge arrester is encompassed by an air box having:

• relative electric permittivity ε r = 1.

2) The cast iron and aluminum, shown in green in Fig. 3.3, cover both the internal air with:

• relative electric permittivity ε r = 1,

• electric conductivity σ = 36.9 MS/m. 3) The porcelain coating part having:

• relative electric permittivity ε r = 1,

SURGE ARRESTER

• electric conductivity σ = 0.

4) The tube of stacked varistor whose electric conductivity is very low:

• relative electric permittivity ε r = 1000,

• electric conductivity is between 10 -2 and 10 -9 S/m.

Computational configurations

The arrester is discretized by a tetrahedron mesh composed of 3,654,533 elements, see Fig. which allows to impose a potential difference between the electrodes. The voltage V s is imposed on one side as shown in Fig. 3.5 and zeros on the other one. Besides, the formulation (E.2) is adopted for this study.

Time study

In the previous sub-section, the geometry, the characteristics of domains, and the mesh are presented. In this part, a time domain study for the surge arrester is presented for both linear and non-linear cases.

Linear case

We begin with the linear case, the purpose of which is to verify that a weak current circulates through the central part of ZnO at a frequency of 50 Hz.

A sinusoidal waveform voltage is applied between the terminals of the conductors (red part) as shown in Fig. 3.5, with a magnitude voltage of V smax = 20 kV for 120 time steps. The total computational time is about 2.4 hours using the CG iterative solver with the Split-Jacobi preconditionner.

Non-linear case

As mentioned before, surge arresters consist of a series or parallel stack of ZnO ceramics. These ceramics have a strongly non-linear V (I) characteristic which is primordial since it is the basis for the use of these ceramics in lightning arresters.

Voltage-current characteristic

The resistive component is divided into three regions due to the conduction mechanisms of the ZnO varistor depending on the temperature denoted by θ and on the magnitude voltage imposed V s as shown in Fig. 3.7. we detail each zone as well as their characteristic:

• Zone 1 has the currents that are highly temperature dependent and the characteristic is non-linear. These currents are low which means that they are mainly capacitive. The resistance decreases as the temperature increases.

• Zone 2 shows temperature independence, while the current is never influenced by temperature. Likewise, the characteristic is non-linear. The current inscreases significantly when the imposed voltage increases.

• Zone 3 represents strong magnitude currents, this means that the V (I) character tends to become ohmic (i.e. linear).

For each zone, there is a mathematical equation describing the behavior of the nonlinear characteristic. Since the zone 2 presents the strong non-linearity, we consider it the following to interpret our simulations in this zone. The latter is described by the following equation

||J||= J 0 exp   - AΦ 3/2 B ||E||   (3.1)
where Φ B is a parameter known by the potential barrier, A and J 0 are constants related to the material, and E is the electric field.

Electro-quasistatic model

Since the ZnO component exhibits a non-linear characteristic, the EQS model should be re-written by introducing a conductivity depending on the electric field. The constitutive law is written in the following form

J = σ(E)E. (3.2)
We recall here's for the EQS model, the electric field E is a curl-free field. An electric scalar potential φ such E = -grad φ -V s grad α is introduced.

In the time domain, the potential formulation to be solved for the non-linear EQS based on the charge conservation law reads

div σ(E)grad φ + ∂ ∂t (εgrad φ) = -div σ(E)V s grad α + ∂ ∂t (εV s grad α) .
(3.3)

Note that the electric conductivity σ defined by the equation (3.3) depends nonlinearly on the electric potential in the regions of ZnO varistor. The weak formulation in the whole domain Ω reads as follows

Ωc σ(E)grad φ • grad φ ′ dΩ c + ∂ ∂t Ω εgrad φ • grad φ ′ dΩ = - Ωc σ(E)V s grad α • grad φ ′ dΩ c - ∂V s ∂t Ω εgrad α • grad φ ′ dΩ.
(3.4)

Consequently, the system to solve is given as

G(φ) + 1 ∆t L X φ i+1 = -V s (t i+1 )G(φ)X α - ∆V s (t i+1 ) ∆t LX α + 1 ∆t LX φ i . (3.5)
where G is the non-linear conductivity matrix and L is the same linear permittivity matrix defined in (1.110).

SURGE ARRESTER

Finally, Newton's scheme is used to solve the non-linear system (3.5), while the CG iterative solver is used to solve the resulting linearized system.

Fitted E-conductivity and V (I)

Based on formula (3.1), we try to fit a non-linear characteristic curve. Taking A = 10 4 , the physical characteristic of a semi-conductor domain Φ B = 0.7 eV, and J 0 is considered as a parameter which allows to adjust the curves. The non-linear conductivity curve is presented in Fig. 3.8a while the voltage-current V (I) curve of the arrester when the magnitude of the voltage varies is presented in Fig. 3.8b. The current is mainly resistive with negligible dielectric losses. We recall here that the voltage-current characteristic is temperature independent.

Current flowing through ZnO

A sinusoidal voltage is applied between the terminals of the surge arrester with a maximum voltage of V smax = 10 kV for 40 time steps. The first two periods of the current crossing the varistor are shown in Fig. 3.9, the non-linear effect can be observed. For the computational time, it takes about 24 minutes for each time step.

Scalar electric potential and electric fields

In this part, the scalar electric potential distribution in the whole domain, the electric field distribution in the non-conductive domains as well as the distribution of the current density in the conductive parts (aluminum and the ZnO varistor) are presented.

Figure 3.9: Total current (A) computed for the first two periods under f = 50 Hz.

Electric potential distribution

In Fig. 3.10, the electric potential in a sectional view (Fig. 3.10a) along the geometry as well as in a 3-D view (Fig. 3.10b) are shown.

It seems that the electric potential varies linearly from the bottom to the top of the arrester, going from 0 V to reach the value maximum of the voltage 10 kV.

Electric field distribution

In Fig. 3.11, the electric field in a sectional view (Fig. 3.11a) along the geometry as well as in a top view (Fig. 3.11b) are shown.

The distribution of the electric field is presented when the magnitude reaches about 120 kV/m in the granular layer in the air tube presented in Fig. 3.1.

Current density distribution

In Fig. 3.12, the current density distribution in a 3-D view is illustrated. We note that the current is concentrated between the highly conductive part and the less conductive part in ZnO. As the section decreases, the current density increases. It is therefore possible to see a channel appear in the center through which almost all the current passes.

In addition, we notice that a large current begins to flow through the surge arrester due to the non-linearity. In other words, the conductivity of a ZnO varistor increases when the voltage at its terminals increases. In fact, we want the surge arrester to become very conductive in the case of over-voltages (for example during a lightning strike).

Synthesis

A complete study on the surge arrester has been carried out. To conclude on the possibilities of modeling surge arresters, surge arresters are used to limit voltage levels. When the voltage becomes too high at the arrester terminals, the varistor ZnO becomes more conductive and the impedance of the arrester decreases. This allows the current to be evacuated to the ground of the electric network. In addition, it would be necessary to introduce a non-linear electric constitutive law on the conductivity of ZnO. In the frequency domain, a study is made with a linear electric constitutive law in order to present the behavior of the surge arrester by computing the global quantity such as the impedance. On the other hand, a spectral approach is essential to take into account the non-linear characteristics.

CHAPTER 3. APPLICATIONS

Industrial applications

Introduction

To validate the Darwin model, two different examples are carried out where the measurements are provided. We are interested in comparing the simulation results

obtained by the Darwin model with the measurements. In the first example, an industrial case represented by a coil designed to characterize the magnetic cores, is studied. For the second one, a complex industrial case composed of two windings insulated and rolled around a toroidal magnetic core, is studied.

Magnetic toroid characterization system

Model presentation

Geometry

The device presented in Fig. 3.14 is composed of two parts: a conductor in golden color of nature OT4228P and a magnetic core in gray color of nature N30-ferrite.

This device is used to characterize the magnetic toroids. A detailed representative scheme is shown in Fig. 3.15. The conductor has two sections, an interior section of 7.07 mm 2 called input surface to enter the current and another exterior of 167.7 mm 2 named output surface to exit it, with a distance of 5 mm between the two sections.

For the magnetic core located in the conductor, has an exterior radius of R = 18 mm, interior radius of r = 11.5 mm, and a height of h = 15 mm.

Material characteristics

The conductive part OT4228P which are very good electric conductors having a conductivity close to that of copper and one magnetic core having a variable permeability as a function of frequency. All the characteristics for all these materials are presented below:

1) The electromagnetic device is encompassed by an air box having a:

• relative magnetic permeability µ r = 1,

• relative electric permittivity ε r = 1.

2) The conductor having a:

• relative magnetic permeability µ r = 1, • relative electric permittivity ε r = 1,

• electric conductivity σ = 15.9 MS/m.

3) The N30-ferrite having a:

• relative magnetic permeability described by the following formula

µ r = µ i 1 + (f /f r ) 2 (3.6)
with µ i the initial permeability equal to 3750 for this material and f r the relaxation frequency equal to 1.8 MHz.

• relative electric permittivity ε r = 5 × 10 4 ,

• electric conductivity σ = 2 MS/m.

Computational configurations

The electromagnetic device is spatially discretized by a tetrahedral mesh composed of 2,866,039 elements. The mesh of the conductor and the magnetic core are presented in Fig. 3.16a and Fig. 3.16b, respectively. The DoFs vector X φ is 477,304 and the DoFs vector X A is 3,341,765. It should be mentioned here that the resistive effects in the conductor are negligible compared to the resistive effects in the toroid when the frequency increases. Consequently, the skin effect should be taken into account in the magnetic core having a thickness of the order of 5 mm at 10 MHz [START_REF] Skutt | Characterization of dimensional effects in ferritecore magnetic devices[END_REF].

INDUSTRIAL APPLICATIONS

Similar to the previous application, the device is enclosed in an air box on which boundary conditions are imposed. Then, a sinusoidal voltage is applied between the terminals of the conductor (gold part) as shown in Fig. 3.15. The frequency interval is [10 : 10 8 ] Hz. The BiCGSTAB is used with a Split-Jacobi preconditioner to solve the system given by invoking the Darwin model ( 2.5) and the computational time for one frequency takes about 1.45 hours.

Measurement and simulation comparison

This model shows not only the eddy current distribution in the magnetic core at high frequencies, but also the displacement currents, which were usually ignored in the quasistatic regimes. The magnetic cores are in fact composed of semi-conductor materials. Due to the small volume of the high frequency magnetic component and the ignorance of the high dielectric constant in ferrite, the dimensional resonance effect was generally overlooked in the field of power electronics. However, with the further increase in switching frequency and power ratings, large magnetic cores can no longer be avoided in some high frequency applications. In Fig. 3.17, a comparison of simulation results and measurement data is shown.

First, it should be mentioned here that the drop in the numerical simulation results at low frequencies is due to the contact resistance which is not taken into account in the modeling.

Moreover, when the frequency increases, the results of the numerical simulation and the measurement show a notable affinity for a large range of frequencies and especially when the capacitive effects become important, in particular, for f > 1 MHz. In addition, the rise in the impedance curve, for f > 20 MHz, does not return to the magnetic toroid, however, it is due to the electromagnetic device which has a greater no-load inductance than the toroid. On the other hand, concerning the phase, it is normal that in the numerical simulation the curve is above since it is less resistive than the other cases. Then, a blockage is noticed at about 90 • , representing a resistive-inductive behavior. This comes down to the fact that the imaginary part µ ′′ of the magnetic permeability µ is not taken into account in the simulation. To validate the results obtained, two models are plotted based on the classic formula

Z 0 = (R 0 + jL 0 ω) + (R core (ω) + jL core (ω)ω) (3.7)
corresponds to the equivalent circuit of the example studied in Fig. 3.14 without taking into account the capacitance, with R 0 and L 0 are the resistance and inductance of the measuring device and R core and L core are that of the magnetic core Figure 3.17: Modulus of impedances and the phase computed with respect to frequencies.

defined as follows

R core = µ 0 µ ′′ (ω) lnN 2 2π ln(R/r)ω, L core = µ 0 µ ′ (ω) lnN 2 2π ln(R/r)
where µ = µ ′ + jµ ′′ with µ ′ and µ ′′ are frequency dependent variables known as the Debye model [START_REF] Hamilton | The small-signal frequency response of ferrites[END_REF] and N is the number of turns which is equal to 1 for this example.

Afterwards, the first model takes into account only the real part µ ′ of µ which shows close results to the simulation and the second one takes into account the two complex parts which shows a drop in the phase such as in the measurement results.

Hence, the drop of the phase observed in Fig. 3.17 for f > 100 kHz is due to the resistive effects which are taken into account by the imaginary part of the magnetic

Toroid core

In this sub-section, we validated the simulation results obtained with Darwin model by considering an industrial transformer. The experimental measurement data, done in a common mode test, are provided until 100 MHz. However, based on the measurement results, it should be mentioned here that the Darwin model is not valid when the frequency is superior to 10 MHz since the radiation effects become non-neglected. As the frequency increases, to handle the skin effect in the periphery of the windings, more elements are needed in the mesh which introduces a scientific challenge at the numerical aspect. For our computation, the total DoFs has reached 40 million elements. The electromagnetic model that we aim to model and study is presented in Fig. 3.18.

Figure 3.18: Toroidal core with two wound windings-side.

Model presentation

Geometry

For the second core example, a single-phase common mode coil consists of two windings wound in the same direction around a toroidal core as shown in Fig. 3.18 and Fig. 3.19a. The windings are made of copper presented in red in Fig. 3.19b, it is circular in shape and each conductor is composed of 10 turns with a section of 0.5027 mm 2 . The enamel component is covering the winding with a thickness of 40 µm. The toroid used for this experiment is a TDK-Epcos toroid reference B64290L0618X830, of material N30 presented in gray in Fig. 3.19 where the material that is in blue is the embedding having a thickness of 0.25 mm.

Material characteristics

The windings are made of copper. The magnetic core has a frequency dependent permeability with different initial assumptions than the previous example. All the materials are presented below:

1) The electromagnetic device is encompassed by an air box having a:

• relative magnetic permeability µ r = 1,

• relative electric permittivity ε r = 1.

INDUSTRIAL APPLICATIONS

2) The conductor having a:

• relative magnetic permeability µ r = 1,

• relative electric permittivity ε r = 1,

• electric conductivity σ = 59.6 MS/m.

3) The insulation covering the winding having a:

• relative magnetic permeability µ r = 1,

• relative electric permittivity ε r = 5.4.

4) The N30-ferrite having a:

• relative magnetic permeability described by the same formula (3.6) with µ i is equal to 4300 for this material and f r is equal to 1.8 MHz,

• relative electric permittivity ε r = 80,

• electric conductivity σ = 2 MS/m.

5) The coating of the magnetic core having a:

• relative magnetic permeability µ r = 1,

• relative electric permittivity ε r = 3.9.

Computational configurations

An electric potential difference φ = V s is imposed between the terminals of the conductor of the windings marked in red as shown in Fig. 3.19b. A single-phase common mode test is considered in this example, which means that φ = V s is imposed on one terminal of the conductor and φ = -V s is on the other as shown in Fig. 3.19a. The frequency interval is [0 : 10 6 ] Hz. To handle the skin effect, the considered mesh is well refined and suitable for all frequencies up to 1 MHz, which features 29,832,477 tetrahedrons including 5,257,323 nodes and 35,093,396 edges.

For each frequency, the computing time using Darwin model takes about 9 days for 40,338,477 DoFs using the BiCGSTAB solver. It is important to point out here that for a system having a huge DoFs, the BiCGSTAB is always preferred than direct solvers because the latter may fail due to the relative need for memory required which is much large than that required for the iterative solvers. The mesh of the winding and the toroidal core are presented in Fig. 3.20.

Figure 3.20: Considered mesh.

Electric potential and electric field

The electric scalar potential along the magnetic core is presented in Fig. 3.21a, which shows a decreasing linear variation of its magnitude. In addition, the electric field distribution in a sectional view is also given in the magnetic core as shown in

Conduction and displacement current densities

When the frequency increases, the total current is divided into two parts: the conduction current in the conductive domains and the displacement current that is located in the non-conductive domains. This means that the displacement current intervenes in the modeling of the electric field in addition to the eddy current. On one hand, the conduction current density is presented in Fig. 3.22a in the windings for a frequency of 500 kHz. High values are observed inside each of the two conductors. Due to the eddy current effects, the current is concentrated in the small layer at the boundary of the conductor. Moreover, the peaks in this case are higher near the middle section mainly due to the magnetic coupling as known the proximity effect.

On the other hand, as shown in Fig. 3.22b, the displacement current density is presented in a sectional view of the magnetic core for a frequency of 500 kHz. The distribution shows that the significant quantity of current is located between the CHAPTER 3. APPLICATIONS turns as well as the magnitude tend to be zero in the contact area with the windings due to the capacitive effects.

Measurement and simulation comparison

In Figure 3.23, the evolution of the impedance as a function of the frequency obtained from the Darwin model as well as the phase are presented.

First, the impedance Z corresponds to the DC resistance of the winding when the frequency tends to zero. Indeed, we have Z = R DC = 6.8 mΩ.

Second, for f < 500 kHz, the results obtained with the Darwin prove a good agreement with the measurement results which present the resistive-inductive phenomena as the MQS regime. MQS regime since the impedance at this point is not linear to the previous points.

The similar behavior can be observed for the phase as a function of the frequency.

Synthesis

In this part, two different industrial applications for Darwin model validation are presented. In the first example, we are interested in showing that the Darwin model shows a good agreement with the measurement results by modeling a simple industrial example. In the second one, the measurement results for an industrial toroidal core are compared with the simulation results by modeling a complex industrial example, which also shows a good agreement, in particular, in the range of intermediate frequencies.

Modeling of an inductance

Introduction

The modeling of the electric devices is currently mainly based on the eddy current problem which allows to obtain a representative model of the behavior in the low frequency domain. With the development of power electronics, devices are subject to high frequency voltage and current stresses that need take into account since this leads in particular to accelerated aging of insulators. Meanwhile, due to the development of the wide bandgap semi-conductors used in the power converters, the frequency of voltage waveforms applied to windings of electric devices becomes increasingly high, especially in the case of pulse-width modulation (PWM).

Consequently, resistive, capacitive, and even inductive effects should be considered simultaneously.

In the literature, the classification or the theoretical limitations of the different quasistatic models are mainly based on the theoretical considerations, in particular, with several assumptions for the computational domain. For example, only the conductive domain is considered, not for the case with multi-domain. The existing reference [START_REF] Koch | Different types of quasistationary formulations for time domain simulations[END_REF] showed a comparison of the results obtained with different quasistatic models only in the time domain while a simple axisymmetric test model represented by a parallel plate capacitor is used. The comparison results were illustrated for the electric field.

The motivation of this work is to investigate the limits of the quasistatic models namely, EQS, MQS, and Darwin models in the general case, in particular, for the computational domain including both conductive and non-conductive domains, which has not been well addressed in the literature. Besides, we are interested in the frequency domain due to the impedances computation. In this example, we compared all the quasistatic models in the frequency domain, by regarding different numerical aspects, namely, the current density, the electric field, the magnetic field, and also the impedance.

Model presentation

In this example, we seek to determine the operating frequency domain of different models by considering an electromagnetic application. Some validation work has been done with the static or/and MQS and EQS models as shown in [START_REF] Steinmetz | Domains of validity of quasistatic and quasistationary field approximations[END_REF][START_REF] Mazauric | Enhancing quasi-static modeling: A claim for electric field computation[END_REF][START_REF] Rapetti | Implications of Galilean electromagnetism in numerical modeling[END_REF]. The objective in this part is to study the influence of the frequency on the distribution of electromagnetics fields, as well as on the impedance obtained by using the different quasistatic models mentioned in the chapter I.

Geometry

As shown in Fig. 3.24, the electromagnetic device composed of a winding (shown in red color) wound around the central column of the magnetic core (shown in blue color) is considered. We have a conductor of section 11.36 mm 2 which takes a rectangular shape, it is composed of 11 turns. The insulation is covering the winding with a thickness of 0.25 mm. The magnetic core is composed of N30-ferrite. This inductance is shown 2-D cutting plane in Fig. 3.25.

Material characteristics

This device is composed by a copper part presented in red, dielectrics part in the grey, and a magnetic core part in blue. The characteristics for all materials are presented below:

1) The electromagnetic device is encompassed by an air box having a:

• relative magnetic permeability µ r = 1,

• relative electric permittivity ε r = 1.

2) The copper having a:

• relative magnetic permeability µ r = 1,

• relative electric permittivity ε r = 1, 3) The insulation covering the winding having a:

• relative magnetic permeability µ r = 1,

• relative electric permittivity ε r = 5.4.

4) The insulation between the magnetic core and the conductor having a:

• relative magnetic permeability µ r = 1,

• relative electric permittivity ε r = 6.4.

5) The N30-ferrite having a:

• relative magnetic permeability µ r = 4300,

• relative electric permittivity ε r = 80,

• electric conductivity σ = 2 MS/m.

Computational configurations

To take into account the skin effects in the winding for high frequencies, the electromagnetic device is spatially discretized by a tetrahedral mesh composed of 14,443,563 elements including 2,466,531 nodes and 16,911,563 edges. The mesh of the winding and of the magnetic core are presented in Fig. 3.26a and Fig. 3.26b, respectively. The Dirichlet condition is imposed on the terminals of the conductor as a potential difference is imposed between the terminals. Then, a sinusoidal voltage is applied

MODELING OF AN INDUCTANCE

between the terminals of the winding marked in red as shown in Fig. 3.24. The electric potential on the below face of the magnetic core is taken to be zero which corresponds to connect to the ground. The frequency interval is [0 : 10 5 ] Hz.

As mentioned before, this application is considered to compare the resulting impedance curves of each model. Then, the computation is done for the EQS, MQS, and Darwin models. For the EQS model, the number of DoFs is 2,466,002 and the computational time for one frequency takes about two hours, while for MQS model it takes about 75 hours with 19,316,839 DoFs. The computational time using the Darwin model takes about 86 hours for 19,373,155 DoFs.

The BiCGSTAB solver is used here with a Split-Jacobi preconditioner to solve the asymmetric system AX = F of the Darwin model (2.5). Consequently, the number of iterations associated for the iterative solver (BiCGSTAB), required to achieve the solution convergence is not the same. For example, for f = 10 kHz and with a fixed stopping criterion, the number of iterations associated for BiCGSTAB is about 14k using the MQS model while it takes about 17k for the Darwin model as shown in Table 3.1. Table 3.1: Comparison between the different models at frequency f = 10 kHz.

Distribution of electric and magnetic fields

At low frequency, when the frequency tends to zero, the distribution of the current density in the winding obtained with all models is similar as observed in Fig. 3.27a, since this is equivalent to solve an electrokinetic problem. In this case, the distribution of the current density is only influenced by the geometry of the winding. When the frequency supplying the winding increases, coupled effects may be observed. To represent these effects in the winding, the distribution of the current density J is presented in Fig. 3.27b-3.27d. At f = 10 kHz, with the EQS model, it is clear that J will never be homogeneous along the winding which means that the skin effect is neglected. On the other hand, the current density with the MQS model is restricted to a small layer at the boundary of the winding, it starts to flow in the outskirts of the conductor, as observed by the eddy current effects. Besides, the coupled for all models, we obtain Z = R DC = 3 mΩ when the frequency tends to zero. When the frequency increases, the eddy current effects appear in the winding, then the value of the resistance of the winding increases. This effect is not taken into account with the EQS model. For f ∈ [0 : 10 3 ] Hz, the MQS and Darwin models give a similar evolution of the impedance. Indeed, the capacitive effect is negligible compared with the inductive and resistive effects. For f > 1 kHz, the influence of the the capacitive effects appears on the evolution of the impedance for both the EQS and Darwin models. A similar behavior can be observed for both models in Fig. 3.30.

Synthesis

At low frequencies, when the impedance exhibits resistive behavior, the EQS model is sufficient. Moreover, when the inductive effects become non-neglected, the MQS and Darwin models can be invoked, but the MQS model is preferred since its computational cost is less expensive than the Darwin model. However, in the intermediate frequencies, in particular, around the resonant frequency, the Darwin model should be adopted since the MQS model cannot handle the coupled capacitive-inductive effects.

Conclusion

In this chapter, a study of an electro-quasistatic case such as the surge arrester is presented in both time and frequency domains. In the following parts, different models have been compared to investigate their limits. The measurement results for an industrial example have been provided to validate the simulation results obtained with the Darwin model, which shows a good agreement in the range of intermediate frequencies, in particular, around the resonant frequency. Additionally, for a complex and multi-domain application, the different quasistatic models have been compared by computing their impedances as well as the different electromagnetic fields in order to investigate their limits.

As EQS and MQS can be seen as the approximation of the Darwin model, we observe numerically that the EQS and Darwin models represent the same behavior when the inductive effects are negligible while the MQS and Darwin models show the same behavior when the capacitive effects are negligible. Furthermore, the EQS model can be used as an indicator to know the frequency when the capacitive effects become significant where the Darwin model should be applied.

Introduction

Due to the development of the wide bandgap semi-conductors used in the power converters, the frequency of voltage and current becomes increasingly high. Then, to study a power converter, the electromagnetic behavior of the printed circuit board (PCB) must be characterized. PCBs can be found in almost all electronic devices and gadgets (cell phones, computers, etc.). Consequently, the PCB industry has rapidly grown due to the great demand.

In Based on a numerical approach such as the finite element method, resistive, capacitive, and inductive effects associated with a PCB should be studied simultaneously when the frequency of excitations increases.

Due to the small size of the tracks, the modeling of a 3-D PCB presents many difficulties in the electromagnetic domain. One of these reasons, the disproportion of the thicknesses of the tracks compared to the other dimensions. The difficulty of dealing with this problem is detailed below.

Numerical issues

The significant contrast existing in the modeling of thin structures such as the PCBs causes difficulties in the generation of meshes. Then, this contrast leads to undesirable elements such as the flat elements. In addition, due to the ill-conditioning of the resulting matrix, the solving becomes more difficult. Besides, another approach can be applied to avoid the meshing of the thin structures.

Generally, the electromagnetic modeling of structures including thin shell model and thin wires is a complex problem in the fields of electrical engineering [START_REF] Rodger | Transient 3D eddy currents in thin sheet conductors[END_REF][START_REF] Krahenbuhl | Thin layers in electrical engineering-example of shell models in analysing eddy-currents by boundary and finite element methods[END_REF][START_REF] Mayergoyz | On calculation of 3-D eddy currents in conducting and magnetic shells[END_REF]. As mentioned before, the geometry of the thin region is characterized by a high ratio between the length and the thickness. Thus, the use of a volume mesh leads to a large number of elements. Furthermore, when the frequency is high the skin depth becomes much smaller than the thickness of the structures. Therefore, the mesh size should be considerably increased to capture the eddy current effects.

Shell elements

The shell element method consists in representing the thin structures by 2-D surfaces, separate two domains meshed with 3-D elements. The use of these elements is imposed at first in the mechanical and thermal problems [START_REF] Surana | Three dimensional curved shell finite elements for heat conduction[END_REF]]. In the 1990s, these elements are introduced into the electromagnetism domain by [START_REF] Brunotte | Line element for efficient computation of the magnetic field created by thin iron plates[END_REF].

The main idea of this method is to avoid the meshing of the thin structures by modeling them in surfaces without thickness. This can be done by developing suitable formulations according to the field continuity conditions through thin structures [START_REF] Rodger | Transient 3D eddy currents in thin sheet conductors[END_REF][START_REF] Krahenbuhl | Thin layers in electrical engineering-example of shell models in analysing eddy-currents by boundary and finite element methods[END_REF][START_REF] Mayergoyz | On calculation of 3-D eddy currents in conducting and magnetic shells[END_REF][START_REF] Guerin | Three dimensional magnetostatic finite elements for gaps and iron shells using magnetic scalar potentials[END_REF]. These formulations can be solved using either the bound-ary element method (BEM) [START_REF] Krahenbuhl | Thin layers in electrical engineering-example of shell models in analysing eddy-currents by boundary and finite element methods[END_REF][START_REF] Pavo | Approximate methods for the calculation of the ect signal of a crack in a plate coated by conducting deposit[END_REF], the finite element method [START_REF] Biro | A finite element formulation for eddy current carrying ferromagnetic thin sheets[END_REF][START_REF] Igarashi | A three dimensional analysis of magnetic fields around a thin magnetic conductive layer using vector potential[END_REF][START_REF] Ren | Degenerated Whitney prism elements-general nodal and edge shell elements for field computation in thin structures[END_REF][START_REF] Geuzaine | Dual formulations for the modeling of thin electromagnetic shells using edge elements[END_REF][START_REF] Choua | Thin crack modeling in ect with combined potential formulations[END_REF] or the integral method such as the PEEC method (Partial Element Equivalent Circuit) [START_REF] Ruehli | Equivalent circuit models for three-dimensional multiconductor systems[END_REF][START_REF] Gope | s)peec: Time-and frequency-domain surface formulation for modeling conductors and dielectrics in combined circuit electromagnetic simulations[END_REF]. Concerning the numerical modeling of the thin regions, a lot of research has been carried out. The applications have been seen for the computation of eddy currents in thin conductive shells [START_REF] Mayergoyz | On calculation of 3-D eddy currents in conducting and magnetic shells[END_REF][START_REF] Biro | A finite element formulation for eddy current carrying ferromagnetic thin sheets[END_REF], for the computation of the magnetostatic field in thin magnetic sheets [START_REF] Biro | A finite element formulation for eddy current carrying ferromagnetic thin sheets[END_REF] and also for the modeling of thin cracks.

Most approaches lead to use a kind of double layer nodal based shell elements [START_REF] Mayergoyz | On calculation of 3-D eddy currents in conducting and magnetic shells[END_REF][START_REF] Guerin | Three dimensional magnetostatic finite elements for gaps and iron shells using magnetic scalar potentials[END_REF], but they are only suitable for interpolating scalar variables. On the other hand, a thin shell eddy current problem is solved using the nodal and the edge version for the discretization of scalar and vector variables, respectively [START_REF] Ren | Degenerated Whitney prism elements-general nodal and edge shell elements for field computation in thin structures[END_REF][START_REF] Koch | Modeling thin conductive sheets using shell elements in magnetoquasistatic field simulations[END_REF][START_REF] Ospina | Electromagnetic field computation in magnetic and conductive thin sheets[END_REF].

In particular, the shell element method provides a more general description of the physical phenomena without any particular hypothesis about the fields inside the thin structure. These elements are built by converting the volume element into a two-dimensional one. In this thesis, a thin shell Darwin problem is solved using the shell element method, applied to both nodes and edges elements. They are derived from the degeneration of Whitney prismatic elements [START_REF] Ren | Degenerated Whitney prism elements-general nodal and edge shell elements for field computation in thin structures[END_REF][START_REF] Koch | Modeling thin conductive sheets using shell elements in magnetoquasistatic field simulations[END_REF] and it can be easily applied to a formulation to solve a thin shell.

Basis functions

The shell elements are obtained by converting the volume element for example the prism or the hexahedron into two-dimensional element as the triangle or rectangle.

Most of the approaches lead to the use of nodal shell elements with a layer of duplicated nodes that allow scalar potentials to be taken into account [START_REF] Guerin | Three dimensional magnetostatic finite elements for gaps and iron shells using magnetic scalar potentials[END_REF][START_REF] Ren | Degenerated Whitney prism elements-general nodal and edge shell elements for field computation in thin structures[END_REF]]. To correctly interpolate the vector potentials, the edge shell elements have been introduced [START_REF] Dular | Mixed finite elements associated with a collection of tetrahedra, hexahedra and prisms[END_REF][START_REF] Ren | Degenerated Whitney prism elements-general nodal and edge shell elements for field computation in thin structures[END_REF].

Moreover, the simultaneous use of the nodal and the edge shell elements allows interpolation of scalar and vector variables in thin structures [START_REF] Ren | Solving 3-D eddy current problem containing thin cracks using dual formulations and shell elements[END_REF][START_REF] Abenius | Thin sheet modeling using shell elements in the finite-element time-domain method[END_REF][START_REF] Koch | Modeling thin conductive sheets using shell elements in magnetoquasistatic field simulations[END_REF].

The shell element is produced by the deformation of a prismatic element as shown in Fig. 4.5. These deformed elements are built by doubling this surface mesh, i.e. both the nodes and the edges of a triangle. In other words, to obtain the shell element, the opposite nodes (1 and 1', 2 and 2', 3 and 3') are brought together to form a surface element where the nodes and the edges are duplicated, but geometrically In the following, the nodal interpolation functions, shown in Fig. 4.6, are defined by

ω 0 i = λ i β, (4.1)
ω 0 i ′ = λ i β ′ , ∀i = 1, 2, 3 (4.2)
where λ i are the reference right triangle functions defined by If W 0 t is the set of nodal functions defined on a tringular mesh then the nodal shell functions can represent the union of the nodal triangular functions defined on the lower and upper side of the prism. We denote in the following by W 0 s the set of the nodal shell functions [START_REF] Ren | Degenerated Whitney prism elements-general nodal and edge shell elements for field computation in thin structures[END_REF] which can be represented by the nodal triangle element W 0 t on the lower and upper surfaces of the prism

             λ 1 = 1 -ξ -η, λ 2 = ξ, λ 3 = η,
W 0 s = (W 0 t β, W 0 t β ′ ). (4.4)
In particular, we have grad β = -n/e and grad β ′ = n/e where n is the unit normal vector of triangles. Then, the gradient of W 0 s reads as follows grad W 0 s = grad s W 0 t β -W 0 t n/e, grad s W 0 t β ′ + W 0 t n/e . w 1 ij = (λ i grad s λ j -λ j grad s λ i )β, ∀i ̸ = j, (4.6)

w 1 i ′ j ′ = (λ i grad s λ j -λ j grad s λ i )β ′ , ∀i ′ ̸ = j ′ , (4.7) w 1 ij ′ = λ i (β ′ gradβ -βgradβ ′ ), ∀i = j ′ . (4.8)
Likewise, the shape functions of the edges defined on a triangle are the terms between the parenthesis of (4.6) and (4.7). Therefore, we denote in the following by W 1 s the set of edge shell functions degenerated from the prism edge element which can be written as follows [START_REF] Ren | Degenerated Whitney prism elements-general nodal and edge shell elements for field computation in thin structures[END_REF]]

W 1 s = (W 1 t β, W 1 t β ′ , W 0 t n/e).
(4.9)

The third term of (4.9) is concluded from (4.8). In fact, since beta β + β ′ = 1 and gradβ = -gradβ ′ = n/e, (4.8) can be simplified as

w 1 ii ′ = λ i n/e.
(4.10)

The rotational of W 1 s reads as follows It should be mentioned here that the prism collapsed into a triangular element while retaining the DoFs associated with the prism as shown in Fig. 4.5-4.7.

curl W 1 s = curl W 1 t β -(n × W 1 t )/e, curl W 1 t β ′ + (n × W 1 t )/e,

Approximation of φ and A

According to (4.4),(4.5) and (4.9),(4.11), the electric scalar potential φ and the magnetic vector potential A, as well as their gradient and their rotational, are approximated by the nodal and edge functions of the shell elements, respectively.

In the shell region, φ is discretized in each shell element by the nodal interpolation function as follows

φ = Nn i,i ′ =1 (βφ i + β ′ φ i ′ )w 0 i . (4.12)
Then, the gradient of φ will be expressed as

grad φ = grad   Nn i,i ′ =1 (βφ i + β ′ φ i ′ ) w 0 i   = Nn i,i ′ =1 grad (βφ i + β ′ φ i ′ )w 0 i + (βφ i + β ′ φ i ′ )grad s w 0 i = Nn i,i ′ =1 n e (φ i ′ -φ i )w 0 i + (βφ i + β ′ φ i ′ )grad s w 0 i (4.13)
where N n is the cardinal of the set of nodes of the triangular mesh.

In addition, it is necessary to compute the circulation over the edges of A, that are represented in the shell elements as

A = Ne i=1 Ne i ′ =1 (βA i + β ′ A i ′ )w 1 ii ′ . (4.14) CHAPTER 4. MODELING OF PCB
Then, the rotationel of A will be expressed as

curl A = curl Ne i=1 Ne i ′ =1 (βA i + β ′ A i ′ ) w 1 ii ′ = Ne i=1 Ne i ′ =1 grad (βA i + β ′ A i ′ ) × w 1 ii ′ + (βA i + β ′ A i ′ )curl s w 1 ii ′ = Ne i=1 Ne i ′ =1 n × w 1 ii ′ (A i ′ -A i )/e + (βA i + β ′ A i ′ )curl s w 1 ii ′ (4.15)
where N e is the cardinal of the set of edges of the triangular mesh.

The approximations defined above, allow a natural continuity of the variables on the interfaces of the thin shell and the various neighboring domains. Then, the A -φ formulation with imposed voltage V s defined in Ω s is written:

The Darwin model in thin structures

           curl s (νcurl s A) + σ ∂A ∂t + grad s φ + ∂ ∂t (εgrad s φ) = -σV s grad s α - ∂V s ∂t (εgrad s α), div s σ ∂A ∂t + grad s φ + ∂V s ∂t (εgrad s φ) = -div s σV s (grad s α) + ∂V s ∂t (εgrad s α) .
(4.16)

The functional spaces of the differential operators defined in Ω s as well as the variational formulations of the system (4.16) are presented in the appendix F. In addition,

APPLICATIONS

it should be mentioned here that by applying shell elements to a variational formulation, the volume integrals concerning the thin shell are transformed into surface integrals where the scalar and vector variables are approximated as (4.12)-(4.15).

These integrals are presented in the appendix G.

Applications

First, we validate initially the implementation of the shell elements by considering two models, one having the conductive track modeled in 3-D and the other modeled in 2-D by using shell elements. The comparison results obtained by the Darwin model are presented in this study. Then, two models of PCBs will be studied in the following parts and the results of the measurements will be provided.

Numerical validation of shell elements

In this part, we present the results of the comparison obtained by different meshes for a simple test case composed of a conductor in L shape enclosing by an air box, as shown in Fig. 4.9, in order to validate the implementation of shell elements.

Material characteristics

The characteristics for the materials are presented below:

123 CHAPTER 4. MODELING OF PCB
1) The conductor is encompassed by an air box having a:

• relative magnetic permeability µ r = 1,

• relative electrical permittivity ε r = 1.

2) The conductor having a:

• relative magnetic permeability µ r = 1,

• relative electrical permittivity ε r = 1,

• electrical conductivity σ = 1 MS/m.

Computational configurations

The track thickness is set at 1.05 mm where it is enlarged by a factor of around 10 with respect to the order of the thicknesses of the tracks. Since the shell elements studied in this thesis are derived from the degeneration of prismatic elements as Table 4.1: Nodal and edge unknowns with the two different meshes.

Since the model is solved by the Darwin formulation, the unknowns of the system are defined on the nodes and the edges as shown in Table 4.1. Then, the DoFs vector X φ has dimension 11,402 for the 3-D mesh with it is equal to 11,063 for the 3-D/2-D mesh and the DoFs vector X A has dimension 71,197 for the 3-D mesh while it is equal to 67,529 for the 3-D/2-D mesh. Note that the difference existing in the DoFs vectors is due to the mesh of the air box which is not exactly the same for the two cases.

The Dirichlet condition is imposed on the terminals of the conductor by imposing a difference of electric potential. A sinusoidal voltage is applied between the terminals of the conductor with V smax = 1 mV at the frequency of f = 100 kHz. The BiCGSTAB is used here with a Split-Jacobi preconditioner to solve the resulting system.

Evolution of the global current

The evolution of the transient currents for the first fives periods is presented with the two cases in Fig. 4.11. The first case, where the mesh is completely meshed by 3-D elements, is represented by the blue curve and the second, where the mesh is made by a coupling of 3-D/2-D elements, is represented by the red curve. The comparison results show that the current is well computed by the shell elements.

Electric and magnetic fields in the track

In this part, we compare the distribution of the current density J and of the magnetic flux density B obtained by the both meshes, as seen in Fig. 4.12 and Fig. 4.13. 1) The electromagnetic device is enclosed by an air box having a:

• relative magnetic permeability µ r = 1,

• relative electrical permittivity ε r = 1.

2) The copper having a:

• relative magnetic permeability µ r = 1,

• relative electrical permittivity ε r = 1,

• electrical conductivity σ = 59.6 MS/m.

3) The epoxy having a:

• relative magnetic permeability µ r = 1,

• relative electrical permittivity ε r = 4.4.

Computational configurations

.2.4 Computation of energies

In this part, we compute the magnetic and electrical energies in order to show that the capacitance is neglected in the studied frequency range. Based on the power balance S, we have S = jωW m + jωW e + P j (4.17)

= V I * (4.18)
where S, V , and I * represent the complex apparent power, the complex voltage and the conjugate of a complex of current, respectively. In addition,

W m = 1 2 L||I|| 2 L 2 (Ω) , W e = 1 2 C||U || 2 L 2 (Ω) , P j = R||I|| 2 L 2 (Ω)
represent the magnetic energy, the electrical energy, and the Joule losses, respectively. Indeed, in Fig. 4.19, the curves of the first two terms of (4.17) are presented.

It is shown that the magnetic energy is dominant and the electrical one is negligible in the frequency interval [1 : 10 6 ] Hz.

Computation of PCB inductance

To compute the PCB inductance, we should first know the wires inductors. Indeed, the global magnetic inductance L g of the system is written as follows

L g = L pcb + L c (4.20)
where L pcb and L c are the inductances of the PCB and the cable, respectively.

The magnetic inductance, is computed in the linear case as function of the magnetic energy W m with

W m = 1 2 Ω B • H dΩ = 1 2 L||I|| 2 L 2 (Ω) . (4.21)
In order to compute the inductance of cables, we consider another geometry of which we reduce or we increase the size of the cables compared to the initial geometry. We keep the same mesh assumptions for the two geometries. In this example, we reduce the size of the cables 3 times. Then, we obtain the first equation

L c = 3L ′ c . (4.22)
For the first example, the length of the cables is set to 75 mm, we had a global inductance of L g = 128 nH. On the other hand, for the second example, the cable length is fixed at 25 mm, hence the new inductance is equal to L ′ g = 52.6 nH. Since the inductance of the PCB will remain constant for both cases, we will have a system of equations with three unknowns Therefore, the inductance of pcb is L pcb = 13.9 nH.

Experimentally, the inductance value of the PCB is 11 nH indicating that the value obtained by the numerical simulation is close. It should be mentioned here that we present the results of the inductance for a single frequency since it is constant in our frequency computation interval.

Electric potential and current density distribution

In Fig. 4.20, the distribution of the electric scalar potential in the whole domain is presented. We notice a linear decrease taking place along the track in a normal way starting from the maximum value of the excitation voltage V smax = 1 V for f = 1 Hz.

Material characteristics

The double-sided model that we aim to study is made up of four different domains.

The 2-D tracks and the 3-D via are made of copper, the epoxy is an insulator that separates the conductive tracks, and the air box. the characteristics of the materials are the same defined for the previous application. CHAPTER 4. MODELING OF PCB

Conclusion

In this chapter, the behavior of a PCB is studied based on the Darwin model in order to capture the coupled resistive, capacitive, and inductive effects. In addition, the thin surfaces are modeled by the shell elements due to the significant contrast existing in the different domains in order to avoid the undesirable elements. Finally, the numerical validations are presented in the section 4.4 as well as two industrial applications on looped PCBs are studied where the measurement results are provided only for the single-sided example and show a good agreement.

Conclusion

With the advent of power electronics, in particular, when the excitation frequency increases, the machines are subjected to high frequency voltage and current stresses that should be taken into account since this leads in particular to accelerated aging of the insulators as well as a breakdown in the dielectrics. Besides, the inclusion of capacitive effects on top of the inductive ones may become necessary.

Since the magneto-quasistatic model does not take into account the capacitive phenomena, the aim of the thesis was to seek and develop tools allowing the modeling of these capacitive and inductive effects simultaneously. The Darwin model which is a surrogate model for full Maxwell system, attracts more and more attention and can be a good choice. This model is able to capture the coupled capacitive-inductive effects in the intermediate frequency row but neglects the radiation effects.

In the first chapter, we returned to the mathematical model formed by the Maxwell relations as well as to the classical assumptions made in the electrical engineering domain. These assumptions make it possible to simplify the complete Maxwell equations into different approximate models, in particular, the electroquasistatic, the magneto-quasistatic, and the Darwin model. The formulation in A and φ potentials are presented. The finite element method is applied in this thesis to solve the different problems in both frequency and time domains. In the case of time discretization, a classical implicit Euler scheme is used. In the last part, numerical validations on an simple test case are presented in order to investigate the importance and the limits of the different quasistatic models, in particular, the Darwin model. formulation requires the use of direct solvers since a saddle point problem is encountered. However, in the case of an industrial problem, the finite element matrix system to be solved generally has a large size, which cannot be handled by direct solvers due to the memory limitation. In this work, a specially designed formulation adopted with iterative solvers is proposed to solve the problem issue from the industrial applications. The solution of the resulting system by an ad-hoc precondition BiCGSTAB iterative solver has proved effective on several industrial cases.

The third chapter concerns the numerical validation of the various quasistatic models as well as some academic and industrial applications are studied. The aim of this chapter is to validate and compare the different models presented in chapter 1. The first application represents a surge arrester classified as a high voltage electro-quasistatic example. The latter is studied in collaboration with EDF for research and development in both time and frequency domains. In the time domain, two studies are performed, a case where the electric constitutive law is linear in order to compute the amplitude of the current which crosses the arrester, and another case where the law is non-linear in order to characterize the voltage-current curve V (I). In the frequency domain, a study is made with a linear electric constitutive law in order to characterize the behavior of the impedance as a function of frequency. In the second part, two industrial applications are studied in order to validate the results and to investigate the limits of the Darwin model. Measurements are provided for both applications and show good agreement with the numerical simulation. The first industrial example represents a magnetic core characterization system composed of a coil and one magnetic core and the second one is composed of two windings wounds in the same direction around a toroidal core.

In the last part, an electromagnetic device represented by an inductance is modeled in order to compare the different quasistatic models namely, the electro-quasistatic, the magneto-quasistatic, and the Darwin model. An investigation is made in order to compare the influence of the distribution of different electromagnetic fields as the function of the frequency where the electromagnetic effects are illustrated. In addition, a comparison of the impedance curves computed by each model is interpreted. Besides, the electro-quasistatic model can be used as an indicator to know the frequency when the capacitive effects become significant where the Darwin model should be applied.

Finally, the last chapter is devoted to the modeling of PCBs in order to compute their impedances. The Darwin formulations are used to describe the behavior of these applications including the illustration of coupled resistive-capacitive-inductive effects. On the other hand, the modeling of thin structures is one of the difficulties encountered in the field of scientific computing. Due to the imbalance between the thicknesses of these structures and the other dimensions, the resulting finite element matrix is ill-conditioned. Otherwise, an extremely fine mesh of the tracks is necessary which is costly in computational time. In this thesis, we resort to the modeling of these structures with surfaces without thickness by using the shell element method, applied to both nodal and edge elements. This method consists in degenerating the prismatic elements into surface triangular elements. This technique as well as the basic spaces and functions in a shell element are studied. For the numerical validation, an academic case with a conductor in L shape is presented. The perspectives raised for this work are manifold. First, as illustrated in this work, the limitation of the direct solvers requires the use of iterative solvers, which can handle more degrees of freedom. Due to the ill-conditioned matrix resulting from the finite element method of the Darwin model, the numerical solution suffers instability for a large range of frequency. The proposed split-Jacobi preconditioner, which is not difficult to be implemented, can ensure a good performance for the iterative solver, but the limitations exist. Besides, efficient and effective preconditioners for high frequencies should be studied in the future work. Second, the modeling of the coupled capacitive-inductive effects in the surge arrester by invoking the Darwin model can be important. As well as, the computation of its impedance provides a comparison with that obtained by the electro-quasistatic model. Since the surge arrester is an electro-quasistatic model, therefore by comparing the electro-quasistatic results with the results obtained by the Darwin model, we can really seek to investigate the limit of the electro-quasistatic model. In addition, a spectral approach is essential to take into account the non-linear characteristics. Third, we are interested in modeling PCBs with multiple ports as well as we could added passive components on the PCB such as inductors, capacitors, and resistors, etc. Moreover, different applications can be considered using shell elements in order to reduce the number of elements and consequently the computational time, for example, the modeling of planar transformers and busbars. In this case, we will be interested on the computation of the inductance variation during switching. Indeed, the magnetic flux distribution will be change due to the modification of the current paths. Finally, the transformers are powered with pulse-width modulation waveform which involves high frequency harmonics on the electromagnetic fields where the modeling of the capacitive effects between the electrical steels can be interesting.

B Static problems

In static regime means that the electromagnetic phenomena are considered independent of time. In this following, the static problems such as electrokinetic, electrostatic, and magnetostatic will be presented.

B.1 Electrokinetic model

In the case of an electrokinetic problem, we seek to compute the distribution of the current density J in a conductor. Then, the Maxwell equations ( 1 The problem comes down to solving the matrix system

GX φ = -V s (t)GX α (B.6)
where G is the same matrix defined in (1.109).

B.2 Electrostatic model

In the case of an electrostatic problem, we seek to compute the distribution of the Let φ h the approximate electric scalar potential of φ defined in the space W 0 (Ω h ) ⊂ H 1 (Ω). The weak formulation (1.86) in W 0 (Ω h ) is written Find φ h ∈ W 0 (Ω h ) such as

(σgrad φ h , grad φ ′ h ) Ωc h + εgrad ∂φ h ∂t , grad φ ′ h Ω h = 0, ∀φ ′ h ∈ W 0 (Ω h ). (C.1)
φ h is written as a linear combination of the nodal functions. Then, the decomposition of φ in its canonical basis is

C.2 Space discretization of the MQS problem

In the same way as EQS, we present the discrete formulations of weak formulations

(1.89) in order to obtain the matrix form. Let φ h the approximate electric scalar potential of φ defined in the space W 0 (Ω c h ) ⊂ H 1 (Ω) and A h the approximate magnetic vector potential of A defined in the space W 1 (Ω h ) ⊂ H(curl , Ω). The weak formulation (1.89) in W 1 (Ω h ) is written

Find (A h , φ h ) ∈ W 1 (Ω h ) × W 0 (Ω c h ) such as ( 1 µ curl A h , curl A ′ h ) Ω h + σ ∂A h ∂t , A ′ h Ωc h + (σgrad φ h , A ′ h ) Ωc h = 0, σ ∂A h ∂t , grad φ ′ h Ωc h + (σgrad φ h , grad φ ′ h ) Ωc h = 0. (C.5)
We introduce the decomposition of A and φ in their canonical basis:

A h = Ne i=1 A i w 1 i , φ h = Nn i=1 φ i w 0 i .
Now, we replace φ ′ h by w 0 j a vector of bases in W 0 (Ω h ) and A ′ h by w 1 j a vector of bases in W 1 (Ω h ), The equation (C.5) reads

D.2 Time discretization of the MQS problem

By applying the implicit Euler scheme (1.118) to (C.8), the matrix system of the MQS model reads:

       M + 1 ∆t S C t 1 ∆t C G             X A i+1 X φ i+1      =        1 ∆t S 0 1 ∆t C 0             X A i X φ i      (D.2)
To obtain a symmetric system, we multiply the second row by ∆t. Then, the MQS system becomes:

      M + 1 ∆t S C t C ∆tG            X A i+1 X φ i+1      =       1 ∆t S 0 C 0            X A i X φ i      . (D.
3)

E Associated excitation source

E.1 EQS model

In order to add the terms introduced for the source in the weak formulations, it is necessary to rewrite the electric field E. We decompose E into two terms, E u is the unknown electric field and E s is the source electric field. The latter allows us to impose a voltage in the expression of E.

Voltage imposition

The field is written

E = E u + E s = -grad φ -V s grad α. (E.1)
Consequently, we will have the matrix system of the EQS model (D.1) written as follows:

G

Current imposition

In the same way, to impose the current with the scalar potential formulation, we use the relation (1.130). This equation is written as follows: (E.4) And the matrix system of the EQS reads as follows:

       G + 1 ∆t L X t α (G t + 1 ∆t L t ) (G + 1 ∆t L)X α X t α (G + 1 ∆t L)X α             X φ i+1 V s i+1      =      0 -I i+1      + 1 ∆t      L X t α L t LX α X t α LX α           X φ i V s i      (E.5)
where G, L, and α h are defined in (1.109), (1.110), and (1.134). In the next parts, the formulation of MQS will be presented to impose a difference of potentials in the conductor.

E.2 MQS model

Similarly, E is written in the following form:

E = E u + E s = - ∂A ∂t -grad φ -V s grad α. (E.6)
Then, the formulation with imposed voltage is written: Consequently, we will have the matrix system of the MQS problem model (D.3) is written as follows:

       
      M + 1 ∆t S C t C ∆tG            X A i+1 X φ i+1      = -V s (t i+1 )      CX α GX α      +       1 ∆t S 0 C 0            X A i X φ i      (E.
H(div s , Ω) = {u ∈ L 2 (Ω s ) | div s u ∈ L 2 (Ω s )}. (F.3)
Consequently, a link exists between the spaces defined above:

∀u

F.3 Discrete variational formulation

In the same way, one quickly presents the discrete variational formulation defined in Ω s . We calculate the discrete formulations of weak forumlations (F.7) in order to make them in matrix form allows to obtain linear systems. Then, let φ h the approximate electric scalar potential of φ defined in the space W 0 (Ω s h ) ⊂ H 1 (Ω s ) and A h the approximate magnetic vector potential of A defined in the space W 1 (Ω s h ) ⊂ H 1 (curl s , Ω s ), where Ω s h is the discretization of the domain Ω s . The weak formulation (F.7) in W 1 (Ω s h ) is written:

G Computation in a shell element

In this part, we present the degeneration of the prismatic elements in special 2-D elements as well as the computations of the finite element terms in a shell element.

G.1 Basic functions in a prism

G.1.1 Nodal functions

The nodal functions of the upper triangle of the prism are defined as follows

             ω 0 1 = λ 1 β, ω 0 2 = λ 2 β, ω 0 3 = λ 3 β. (G.1)
In the lower triangle of the prism are defined as follows

             ω 0 1 ′ = λ 1 β ′ , ω 0 2 ′ = λ 2 β ′ , ω 0 3 ′ = λ 3 β ′ .
(G.2)

In the following, we compute the gradient of these functions such that

G.1.2 Edge functions

The edge functions of the upper triangle of the prism are defined as follows In the lower triangle, they are defined as follows

            
             w 1 1 ′ 1 ′ = (λ 1 grad s λ 2 -λ 2 grad s λ 1 )β ′ , w 1 2 ′ 2 ′
= (λ 1 grad s λ 3 -λ 3 grad s λ 1 )β ′ , w 1 3 ′ 3 ′ = (λ 2 grad s λ 3 -λ 3 grad s λ 2 )β ′ .

(G.5) For the edges perpendicular to the surfaces, they are defined as follows In the following, an edge function and its rotational on the edge p 1 p 2 for example is written

             w 1 1 ′ 1 ′ = λ 1 /e,
w 1 12 =      (λ 1 + λ 2 )β λ 2 β 0      curl w 1 12 =      λ 2 λ 3 -1 2β      . (G.7)

G.2 Computation of the transformation matrix

We consider a prism defined in the reference axis (1,2,3) and perpendicular along We assume that the prism points have the following coordinates: p 1 (0; 0; -0.5), p 2 (1; 0; -0.5), p 3 (0; 1; -0.5), p 1 ′ (0; 0; 0.5), p 2 ′ (1; 0; 0.5), p 3 ′ (0; 1; 0.5). Then, we have the following system

G.3 Gradient in the real element

By denoting ω0 a nodal function in the real coordinate system, one writes grad ω0 i = (J t ) -1 grad ω 0 i , ∀i = 1, 2, 3. (G.10)

For example, for the gradient of ω0 1 is written as follows (G.12)

grad ω0 1 = (J t ) -1 grad ω 0 1 = 1 |J ∆ |   |J ∆ |(J t
Similarly, in the lower triangle, we have

(J t ) -1 grad ω0 i ′ =   (J t ) -1 grad ∆ ω 0 1 β ′ (J t ) -1 grad ∆ ω 0 2 β ′ (J t ) -1 grad ∆ ω 0 3 β ′ λ 1 /e λ 2 /

F

  characteristics . . . . . . . . . . . . . . . . 4.4.1.2 Computational configurations . . . . . . . . . . . . . 4.4.1.3 Evolution of the global current . . . . . . . . . . . . 4.4.1.4 Electric and magnetic fields in the track . . . . . . . 4.4.1.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . 4.4.2 PCBs loop modeling . . . . . . . . . . . . . . . . . . . . . . . 4.4.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . 4.4.2.2 Single-sided PCB . . . . . . . . . . . . . . . . . . . . 4.4.2.3 Double-sided PCB . . . . . . . . . . . . . . . . . . . 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusion Appendices A Ostrogradski and Green theorems A.1 Ostrogradski's theorem . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Green formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B Static problems B.1 Electrokinetic model . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Electrostatic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Magnetostatic model . . . . . . . . . . . . . . . . . . . . . . . . . . . C Discrete variational formulations in space C.1 Space discretization of the EQS problem . . . . . . . . . . . . . . . . C.2 Space discretization of the MQS problem . . . . . . . . . . . . . . . . D Discrete variational formulations in time D.1 Time discretization of the EQS problem . . . . . . . . . . . . . . . . D.2 Time discretization of the MQS problem . . . . . . . . . . . . . . . . E Associated excitation source E.1 EQS model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.2 MQS model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Functional spaces and variational formulations in the shell structures F.1 Functional spaces of differential operators . . . . . . . . . . . . . . . . F.2 Continuous variational formulation . . . . . . . . . . . . . . . . . . . CONTENTS F.3 Discrete variational formulation . . . . . . . . . . . . . . . . . . . . . G Computation in a shell element G.1 Basic functions in a prism . . . . . . . . . . . . . . . . . . . . . . . . G.1.1 Nodal functions . . . . . . . . . . . . . . . . . . . . . . . . . . G.1.2 Edge functions . . . . . . . . . . . . . . . . . . . . . . . . . . G.2 Computation of the transformation matrix . . . . . . . . . . . . . . . G.3 Gradient in the real element . . . . . . . . . . . . . . . . . . . . . . .List of Figures1.1 Studied domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Domain study of EQS problem. . . . . . . . . . . . . . . . . . . . . . 1.3 Domain study of MQS problem. . . . . . . . . . . . . . . . . . . . . . 1.4 Example of a triangular mesh made up of several domains. . . . . . . 1.5 Discretization of the time interval [0, T ] in i constant time steps . . . 1.6 Example of γ path . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 The geometry of the electric cable. . . . . . . . . . . . . . . . . . . . 1.8 Considered mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Magnitude of the conduction current density J (A.m -2 ) for different frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .(a) f = 10 Hz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) f = 10 kHz . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) f = 100 MHz . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.10 Magnitude of the magnetic flux density B (T) for different frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) f = 10 Hz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) f = 10 kHz . . . . . . . . . . . . . . . . . . . . . . . . . . . . (c) f = 100 MHz . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.11 Magnitude of the displacement current density ∂D ∂t (A.m -2 ) for different frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) f = 10 Hz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) f = 10 kHz . . . . . . . . . . . . . . . . . . . . . . . . . . . . LIST OF FIGURES (c) f = 100 MHz . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.12 Modulus of impedance and the phase computed with respect to frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.13 Electrical power (W) and the relative error (%) are computed as the function of the frequency. . . . . . . . . . . . . . . . . . . . . . . . . . (a) Electrical power . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Relative error . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.14 Diagram specifying the different electromagnetic models. . . . . . . . 2.1 (a) The geometry of the capacitor; (b) the mesh of the geometry. . . (a) Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Variation of the currents w.r.t. the frequency using BiCGSTAB iterative solver (red) and Mumps direct solver (blue) for all formulations (2.5), (2.7), and (2.10). . . . . . . . . . . . . . . . . . . . . . . . . . . (a) In [0 : 10 9 ] Hz frequency interval. . . . . . . . . . . . . . . . . (b) In [10 7 : 10 9 ] Hz frequency interval. . . . . . . . . . . . . . . . 2.3 Computational time of the solving for differents formulations/solvers w.r.t. frequency (BiCGSTAB (red); Mumps (blue)). . . . . . . . . . . 2.4 Number of iterations for differents formulations w.r.t. frequency. . . . 2.5 Relative residue computed by the BiCGSTAB-iterative solver with a Split-Jacobi preconditioner. . . . . . . . . . . . . . . . . . . . . . . . (a) For f = 100 Hz. . . . . . . . . . . . . . . . . . . . . . . . . . . (b) For f = 10 MHz. . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Evolution of the current (A) obtained for the first two periods by the M1, M2, and M3 formulations. . . . . . . . . . . . . . . . . . . . . . . 3.1 The upper part of the geometry. . . . . . . . . . . . . . . . . . . . . . 3.2 Cross-sectional view of the surge arrester. . . . . . . . . . . . . . . . . 3.3 Schematic sectional view. . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Mesh of the studied device. . . . . . . . . . . . . . . . . . . . . . . . 3.5 Section where the potential V s is imposed. . . . . . . . . . . . . . . . 3.6 Evolution of the current (A) flowing through the ZnO varistor. . . . . 3.7 Voltage-current characteristic (source: EDF internal report). . . . . . 3.8 Non-linear characteristics in the surge arrester. . . . . . . . . . . . . . (a) E-conductivity in the ZnO varistor . . . . . . . . . . . . . . . (b) Voltage-Current . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Total current (A) computed for the first two periods under f = 50 Hz. LIST OF FIGURES 3.10 Distribution of the electric potential (V) at t = 25 ms. . . . . . . . . (a) Cross-sectional view . . . . . . . . . . . . . . . . . . . . . . . (b) 3-D view . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Distribution of the electric field (V/m) at t = 25 ms. . . . . . . . . . (a) Cross-sectional view . . . . . . . . . . . . . . . . . . . . . . . (b) Top view . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Distribution of the current density J (A.m -2 ) at t = 25 ms. . . . . . . 3.13 Modulus of impedance and the phase computed with respect to frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Geometry of the electromagnetic device. . . . . . . . . . . . . . . . . 3.15 Schematic representation of the electromagnetic device. . . . . . . . . 3.16 Considered mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Conductor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Magnetic core . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.17 Modulus of impedances and the phase computed with respect to frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Toroidal core with two wound windings-side. . . . . . . . . . . . . . . 3.19 3-D model of toroidal core. . . . . . . . . . . . . . . . . . . . . . . . . (a) Geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Detail of winding. . . . . . . . . . . . . . . . . . . . . . . . . . 3.20 Considered mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.21 Distribution of the electric scalar potential and the electric field in the magnetic core for f = 500 kHz. . . . . . . . . . . . . . . . . . . . (a) Electric scalar potential (V) . . . . . . . . . . . . . . . . . . . (b) Electric field E (V/m) . . . . . . . . . . . . . . . . . . . . . . 3.22 Distribution of the conduction current density and the displacement current density in the windings and the magnetic core, respectively. . (a) Conduction current density J (A/m 2 ) . . . . . . . . . . . . . . of impedance and the phase computed with respect to frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.24 Geometry of the electromagnetic device. . . . . . . . . . . . . . . . . 3.25 Schematic representation of the electromagnetic device. . . . . . . . . 3.26 Considered mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Winding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Magnetic core . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.27 Distribution of the current density J (A.m -2 ) in the winding. . . . . . LIST OF FIGURES (a) All models with f = 1 Hz . . . . . . . . . . . . . . . . . . . . (b) EQS model with f = 10 kHz . . . . . . . . . . . . . . . . . . . (c) MQS model with f = 10 kHz . . . . . . . . . . . . . . . . . . (d) Darwin model with f = 10 kHz . . . . . . . . . . . . . . . . . 3.28 Distribution of the magnetic flux density B (T) in the magnetic core at f = 1 Hz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.29 Magnitude of the electric field E (V/m) between the gaps of the winding at f = 10 kHz. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.30 Modulus of impedance computed with respect to the frequency. . . . 4.1 Printed circuit board. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Example of PCB model. . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Considered meshes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Coarse mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Fine mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Convergence with the different meshes. . . . . . . . . . . . . . . . . 4.5 Right triangular prism element (left) and the corresponding shell element (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Nodal shell element. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Edge shell element. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Studied domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Studied domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Considered meshes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) 3-D mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) 2-D mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Evolution of the current (mA) when the frequency of the excitation is f = 100 kHz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Distribution of the current density J (mA.m -2 ) in the track at t = 35 µs.126 (a) 3-D mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) 2-D mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Distribution of the magnetic flux density B (mT) in the track at t = 35 µs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) 3-D mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) 2-D mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Real device. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 PCB model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 The imposed voltage sections. . . . . . . . . . . . . . . . . . . . . . . 4.17 Considered mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.18 PCB connexion with two electric wires. . . . . . . . . . . . . . . . . . (a) Far view . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) Near view . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.19 Magnetic and electrical energies (J) vs frequency. . . . . . . . . . . . 4.20 Magnitude of the scalar electric potential (V) in the whole domain. . 4.21 Distribution of the current density J (A.m -2 ) in the track for different frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) For f = 1 Hz . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) For f = 100 kHz . . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Modulus of impedance computed with respect to the frequency. . . . 4.23 Voltage source (mV) (square wave). . . . . . . . . . . . . . . . . . . . 4.24 Evolution of the current (mA) when the frequency of the excitation is f = 100 kHz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.25 Real device. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.26 Considered model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) PCB model . . . . . . . . . . . . . . . . . . . . . . . . . . . . (b) 3-D conductive via . . . . . . . . . . . . . . . . . . . . . . . . 4.27 The imposed voltage sections. . . . . . . . . . . . . . . . . . . . . . . 4.28 Considered mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Mesh for whole domain . . . . . . . . . . . . . . . . . . . . . . (b) Mesh in the conductive via . . . . . . . . . . . . . . . . . . . . 4.29 Magnitude of the scalar electric potential (V) in the whole domain. . 4.30 Distribution of the current density J (A.m -2 ) in the tracks and the conductive via. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.31 Modulus of impedance and the phase computed with respect to frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . G.1 Transformation of the real element to the reference element. . . . . . List of Tables 1.1 Electromagnetic models classified according to the physical phenomena taken into account. . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Different cases for the studied formulations. . . . . . . . . . . . . . . 2.2 DoFs for the different cases. . . . . . . . . . . . . . . . . . . . . . . . 2.3 The performance of iterative solver (BiCGSTAB-M3) with a Split-Jacobi preconditioner and direct solver (Mumps-M2) versus the number of elements for different frequencies. . . . . . . . . . . . . . . . . (a) For f = 10 kHz. . . . . . . . . . . . . . . . . . . . . . . . . . . (b) For f = 100 kHz. . . . . . . . . . . . . . . . . . . . . . . . . . (c) For f = 1 MHz. . . . . . . . . . . . . . . . . . . . . . . . . . . (d) For f = 10 MHz. . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 The performance of iterative solver (BiCGSTAB-M1)-reference and (CG-M3) with a Split-Jacobi preconditioner and direct solver (Mumps-M2) versus the number of elements for f = 10 MHz. . . . . . . . . . . 3.1 Comparison between the different models at frequency f = 10 kHz. . 4.1 Nodal and edge unknowns with the two different meshes. . . . . . . .

CHAPTER

  second one entitled "Comparison of different Formulations/Solvers of the Darwin model" is devoted on the methods of solving of the various formulations of Darwin. A state of the art on the methods of solving of linear partial differential equations is detailed in the first part. In the second part and by considering different formulations of the Darwin model, a detailed solving comparison is presented in both frequency and time domains and by using both direct and iterative solvers.The third one entitled "Applications" is devoted on academic and industrial applications. Four applications are carried out in this chapter. The first application, a surge arrester exhibiting a non-linear conductivity law, is studied within the framework of collaboration with EDF for research and development. In the second part, two industrial applications are presented. The first example represents a magnetic core characterization system composed of a coil and one magnetic core while the second example is composed of two windings wound in the same direction around a toroidal core. These two applications are studied in order to investigate the limit and to validate the results of the Darwin model. In addition, the comparison of the simulation results with the measurement results is provided. Finally, in the last part, an academic electromagnetic device represented by an inductance model is studied in order to compare the impedance curves computed by the different quasistatic models, namely EQS, MQS, and Darwin models.The last chapter entitled "Modeling of PCB" is devoted on the study of PCBs in order to characterize the behavior of their impedances. However, due to the disproportion between the conductive tracks and the other dimensions existing in a 3-D model, many difficulties are presented in the numerical simulation. Besides, the modeling of thin structures requires the optimization of the mesh with the shell elements. In this INTRODUCTION chapter, the Darwin model is presented with the shell elements in order to study the complete behavior of the PCBs when the frequency increases. Furthermore, different geometries of PCBs are studied as well as the measurements are provided to validate the numerical simulation results. The presented Ph.D. thesis has been done in the L2EP laboratory (Laboratory of Electrical Engineering and Power electronics of Lille), in the Numerical Tools and Methods team which is co-funded by the «Hauts de France» region within the framework of FEDER (Fonds Européen de Développement Régional) and by «EDF R&D». The work in the thesis is focused on the implementation of the Darwin model using the finite element method in order to study the resistive-inductive-capacitive effects. From a numerical point of view, the code_Carmel is a 3-D finite element computation code developed within LAMEL laboratory, a common laboratory between L2EP laboratory and EDF since 2006. The code is dedicated to model low frequency electromagnetic devices. This amounts to solving the Maxwell equations within the framework of the quasi-stationary regime approximation and taking into account the MQS problem. The code_Carmel is based on the solving of the formulations in potentials (vectors and scalars) using the finite element method. Therefore, one part of the thesis is devoted to implementing the tools in the code_Carmel in order to quantify the electric field in dielectrics such as the electrostatic, the EQS as well as the Darwin model which is adapted to simultaneously compute the electric and magnetic fields at the intermediate frequencies. 1 Modeling of electromagnetic problems Electromagnetic systems are governed by Maxwell's equations. The goal of this first part is to present the formulations with magnetic vector potential and electric scalar potential in three dimensions using the finite element method. Numerical validations are presented in order to compare the different quasistatic models and to determine the need to incorporate the modeling of the coupled capacitive-inductive effects, in particular, the Darwin model. Numerical validation . . . . . . . . . . . . . . . . . . . . . 51 1.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

Fig. 1 .

 1 Fig. 1.1 presents the study domain Ω with the boundary Γ having two complementary parts such as Γ B and Γ H where Γ B ∪ Γ H = Γ and Γ B ∩ Γ H = ∅. The boundary conditions of magnetic quantities are such that

CHAPTER 1 .

 1 MODELING OF ELECTROMAGNETIC PROBLEMS electric called electro-quasistatic (EQS) model, one is magnetic called magnetoquasistatic (MQS) model, and the last one is electric and magnetic at the same time called Darwin model. Let us introduce the velocity v of the system with modulus |v|= L/T , where L and T denote the units of space and time, respectively, and we compare it with the light celerity in the continuous domain c = 1/ √ εµ with c = 3 × 10 8 m.s -1 . The approximation of quasistatic regimes corresponds to |v|≪ c but this condition is not sufficient to choose between the EQS, MQS, and Darwin models [Rapetti and Rousseaux, 2011]. Let us introduce E = E ′ E and B = B ′ B, where E ′ , B ′ are reference quantities and E, B are non-dimensional quantities of order 1. It should be mentioned here that in dimensional analysis, the spatial and time differentiation are equivalent to

Figure 1 . 2 :

 12 Figure 1.2: Domain study of EQS problem.

  ics is a means of studying magnetic effects in an unsteady state. The phenomena considered in this model depend on time. As shown in Fig. 1.3, the MQS problem to be solved may include an air box in addition to different regions, a conductive magnetic region represented by Ω c domain, and a non-conductive magnetic region such as the iron.

Figure 1 . 3 :

 13 Figure 1.3: Domain study of MQS problem.

Figure 1 . 4 :

 14 Figure 1.4: Example of a triangular mesh made up of several domains.

CHAPTER 1 .

 1 MODELING OF ELECTROMAGNETIC PROBLEMS [0, T ] the study interval which is divided into i time step as shown in Fig. 1.5.

Figure 1 . 5 :

 15 Figure 1.5: Discretization of the time interval [0, T ] in i constant time steps

Figure 1

 1 Figure 1.6: Example of γ path

  Most of the work presented in this thesis is done with code_Carmel. The problems considered by the code are focused on the studies of inductive or inductive-resistive coupled effects via the study of the magnetostatic and the MQS problem, respectively. Hence the objective of this part is to validate the systems implemented in the code_Carmel such as the EQS problem and the Darwin model.In this part, we are interested in the validation of the numerical results of the Darwin model, as well as the determination of the limit of the different models studied above by comparing the impedance curves computed by each model. Thus, this academic application dedicates in the validation of the parameters k 1 , k 2 , and k 3 defined in section 1.2.1 and shows the need to use the Darwin model when the frequency increases.

  Figure 1.8: Considered mesh.

  Figure 1.10: Magnitude of the magnetic flux density B (T) for different frequencies.

  Figure 1.11: Magnitude of the displacement current density ∂D ∂t (A.m -2 ) for different frequencies.

Figure 1

 1 Figure 1.13: Electrical power (W) and the relative error (%) are computed as the function of the frequency.

  an academic example is considered to compare the different models as well as a diagram that specifies the different electromagnetic domains including the Darwin model is presented. 2 Comparison of different Formulations/Solvers of Darwin model In this chapter, a state of the art on the methods to solve linear partial differential equations is introduced. Two families of linear solvers are presented in the first part: direct methods and iterative methods. In the second part, a complete and detailed comparison considering different formulations of Darwin model and different solvers is presented.

Contents 2 .

 2 1 Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 2.2 Different formulations/cases . . . . . . . . . . . . . . . . . 63 2.3 Comparison of formulations with different solvers in the frequency domain . . . . . . . . . . . . . . . . . . . . . 69 2.4 Performance of M3 formulation . . . . . . . . . . . . . . . 72 2.5 Comparison of solvers versus number of elements . . . . 75 2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 CHAPTER 2. COMPARISON OF DIFFERENT FORMULATIONS/SOLVERS OF DARWIN MODEL

  of linear equations and it is the most used among the different Krylov methods. It was pioneered by Magnus Hestenes, Eduard Stiefel, and Cornelius Lanczos in the 1950s. It was named as one of the top 10 algorithms of the 20th century. It is considered first and foremost as a direct method. John Reid's work in the early 1970s brought renewed attention to the algorithm, and since then it has been a subject of intense research. In 1971, it was rediscovered by Reid as an iterative method.

  summarize, we consider three different formulations for the Darwin model derived from the Maxwell equations. The first formulation with two unknowns (A, φ) (edge unknowns + nodal unknowns) noted M1 where the resultant FE matrix is ungauged and asymmetric. The second formulation noted M2 having three unknowns (A, φ, p) (edge unknowns + 2 × nodal unknowns) with a symmetric resulting matrix and gauged with the Coulomb technique. Finally, by using the method of double Lagrange multipliers, we adapt the M2 formulation to M3 formulation in order to CHAPTER 2. COMPARISON OF DIFFERENT FORMULATIONS/SOLVERS OF DARWIN MODEL use an iterative solver in the case of large size problems. This formulation have four unknowns (A, φ, p 1 , p 2 ) (edge unknowns + 3 × nodal unknowns).

Figure 2 .

 2 Figure 2.1: (a) The geometry of the capacitor; (b) the mesh of the geometry.

Table 2 . 2 :

 22 DoFs for the different cases. is proportional to the input voltage frequency. Indeed, the current increase linearly as a function of the frequency as shown in Fig. 2.2 (black curve). As well, the results obtained by different cases are very close in the interval [1 : 10 7 ] Hz, showing

  (a) In [0 : 10 9 ] Hz frequency interval. (b) In [10 7 : 10 9 ] Hz frequency interval.

Figure 2 . 2 :

 22 Figure 2.2: Variation of the currents w.r.t. the frequency using BiCGSTAB iterative solver (red) and Mumps direct solver (blue) for all formulations (2.5), (2.7), and (2.10).

Figure 2 . 3 :

 23 Figure 2.3: Computational time of the solving for differents formulations/solvers w.r.t. frequency (BiCGSTAB (red); Mumps (blue)).

Figure 2 . 4 :

 24 Figure 2.4: Number of iterations for differents formulations w.r.t. frequency.

Figure 2 . 5 :

 25 Figure 2.5: Relative residue computed by the BiCGSTAB-iterative solver with a

  (a) For f = 10 kHz. (b) For f = 100 kHz. (c) For f = 1 MHz. (d) For f = 10 MHz.

CHAPTER 2 .

 2 COMPARISON OF DIFFERENT FORMULATIONS/SOLVERS OF DARWIN MODEL

Figure 2 .

 2 Figure 2.6: Evolution of the current (A) obtained for the first two periods by the M1, M2, and M3 formulations.

3 Applications

 3 In this chapter, different academic and industrial applications are modeled and studied in order to compare the different quasistatic models, namely EQS, MQS, and the Darwin model. The first application, a surge arrester exhibiting a non-linear conductivity law, is studied within the framework of collaboration with EDF for research and development. To validate the results and investigate the limit of the Darwin model, two industrial cases are modelled and the results are compared with the measurements. The first industrial example represents a magnetic core characterization system composed of a coil and a magnetic core while the second one is composed of two windings wound around a toroidal core. In the last part, an academic electromagnetic device represented by an inductance model is studied in order to investigate the influence of the frequency on the electromagnetic field distributions and to compare the impedance curves computed by each model. Modeling of an inductance . . . . . . . . . . . . . . . . . . 105 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 3.1 Introduction Before studying a Darwin model, we start our analysis with a complete study on an EQS example with only the capacitive effects. Accordingly, in the first application, a surge arrester application based on the EQS model taking into account the non-linear electric constitutive law is presented. For the second part, two applications are presented to validate the results obtained from the numerical simulation of the Darwin model with measurement results. It is followed in the next part by an electromagnetic device which is studied in order to compare the resulting electromagnetic effects as a function of the frequency. Both the local field distributions and the global impedances are presented to show the necessity of the Darwin model when the frequency increases.

  other. The number and dimensions of these pads vary depending on the nominal voltage of the surge arrester as shown in Fig. 3.1 where only the upper part of the geometry is presented. Their diameter varies between 25 and 80 mm while their height is from 15 to 40 mm. The total height of the surge arrester is 80 cm while it is equal to 70 cm for the ZnO varistor, as shown in Fig. 3.2.

Figure 3

 3 Figure 3.1: The upper part of the geometry.

Figure 3 . 2 :

 32 Figure 3.2: Cross-sectional view of the surge arrester.

Figure 3 . 3 :

 33 Figure 3.3: Schematic sectional view.

  3.4. Since the model is solved by the EQS model, the unknowns of the system are defined only on the nodes and the DoFs is 462,207.

Figure 3 . 4 :

 34 Figure 3.4: Mesh of the studied device.

Figure 3 . 5 :

 35 Figure 3.5: Section where the potential V s is imposed.

Figure 3

 3 Figure 3.6: Evolution of the current (A) flowing through the ZnO varistor.

Fig. 3

 3 Fig. 3.6 shows the current flowing in the ZnO varistor structure computed with the

Figure 3

 3 Figure 3.7: Voltage-current characteristic (source: EDF internal report).

Fig. 3 .

 3 Fig. 3.7 shows the voltage-current curve of the ceramic component. In the following,

  Figure 3.8: Non-linear characteristics in the surge arrester.

Figure 3

 3 Figure 3.10: Distribution of the electric potential (V) at t = 25 ms.

  Figure 3.11: Distribution of the electric field (V/m) at t = 25 ms.

Figure 3

 3 Figure 3.13: Modulus of impedance and the phase computed with respect to frequencies.

Figure 3 .

 3 Figure 3.14: Geometry of the electromagnetic device.

Figure 3 .

 3 Figure 3.15: Schematic representation of the electromagnetic device.

Figure 3 .

 3 Figure 3.16: Considered mesh.

  Fig. 3.19.

Figure 3 .

 3 Figure 3.19: 3-D model of toroidal core.

Fig. 3 .

 3 Fig. 3.21b. It can be observed that the electric field is located in the periphery of the core due to the eddy current effects that occurs and reaches a maximum value of 12 V/m.

  Figure 3.21: Distribution of the electric scalar potential and the electric field in the magnetic core for f = 500 kHz.

2 )

 2 Figure 3.22: Distribution of the conduction current density and the displacement current density in the windings and the magnetic core, respectively.

Figure 3 .

 3 Figure 3.23: Modulus of impedance and the phase computed with respect to frequencies.

  Figure 3.24: Geometry of the electromagnetic device.

Figure 3 .

 3 Figure 3.25: Schematic representation of the electromagnetic device.

Figure 3 .

 3 Figure 3.26: Considered mesh.

  Figure 3.27: Distribution of the current density J (A.m -2 ) in the winding.

  3.4. MODELING OF AN INDUCTANCEinductive-capacitive effects can be observed with the Darwin model. Therefore, due to the capacitive effects, the current density distribution obtained with the Darwin model is not the same as that of MQS.

Figure 3 .

 3 Figure 3.28: Distribution of the magnetic flux density B (T) in the magnetic core at f = 1 Hz.

Figure 3 .

 3 Figure 3.29: Magnitude of the electric field E (V/m) between the gaps of the winding at f = 10 kHz.

  Fig. 4.1, an example has been shown: chips, capacitors, coil, hardware part. A PCB is made up of thin strips of a conductive material such as copper, which have been etched from a layer attached to a flat insulating sheet, and onto which integrated circuits and other components are attached.

Figure 4 . 1 :

 41 Figure 4.1: Printed circuit board.

CHAPTER 4 .

 4 MODELING OF PCBFor these reasons, an example of a PCB proposed in Fig.4.2 is studied in order to highlight some difficulties encountered during the numerical simulation.

Figure 4 . 2 :

 42 Figure 4.2: Example of PCB model.

Figure 4

 4 Figure 4.3: Considered meshes.

Figure 4

 4 Figure 4.4: Convergence with the different meshes.

Figure 4

 4 Figure 4.5: Right triangular prism element (left) and the corresponding shell element (right).

  β = 1/2 -ζ, β ′ = 1/2 + ζ are the linear interpolation functions along the axis ζ.

Figure 4

 4 Figure 4.7: Edge shell element.

  Similar for the Darwin model defined in Ω in (1.61), nothing changes in Ω s , all terms are kept. The same Darwin system is defined in Ω s in addition to the overall system defined in Ω. As represented in Fig.4.8, the problem to be solved includes the Ω ns domain containing the 3-D meshed regions and the Ω s domain containing the 2-D meshed shell surfaces.

Figure 4 . 8 :

 48 Figure 4.8: Studied domain.

Figure 4 . 9 :

 49 Figure 4.9: Studied domain.

Figure 4 .

 4 Figure 4.10: Considered meshes.

Fig. 4 .

 4 Fig. 4.12a shows the distribution of J in the track in 2-D and Fig. 4.12b shows the distribution of J in the 3-D one. A good distribution along the track in 2-D is obtained as well as at the magnitude of J, the results are very close comparing with the 3-D case. Likewise, regarding to the distribution of the magnetic flux density,

Figure 4

 4 Figure 4.12: Distribution of the current density J (mA.m -2 ) in the track at t = 35 µs.

Figure 4

 4 Figure 4.13: Distribution of the magnetic flux density B (mT) in the track at t = 35 µs.

  In this section, the shell element method is presented as well as the associated basis functions for both nodal and edge elements. The implementation of this method is validated by considering a L shape conductive track and by comparing both the global and local results. In the next sections, two applications of PCBs are studied by incorporating the shell elements. inductance has two main origins: the internal interconnections in devices packages and the external interconnections, for example, on PCB tracks. The advantage of the modeling of the looped PCBs is that there is only one component, the copper. Furthermore, the passive and transient components are not taken into account in the modeling. Our interest in modeling a PCB is to compute its impedances.In this section, we study different examples of PCBs having a simple structure like loop PCBs in order to capture the existing electromagnetic effects when the excitation frequency increases. Two different PCB designs of the commutation loop are studied in the following, a single-sided PCB and a double-sided PCB. For the first example, the global and local results obtained by the numerical simulation are presented as well as the comparison with measurements is provided. For the second one, the interest is to study the coupling of 2-D and 3-D domains in conductivity and to describe the capacitive-inductive effects at high frequencies. In the following, the numerical simulation results for the two studied PCBs are obtained by applying the Darwin model in the frequency domain.4.4.2.2 Single-sided PCBFor the first application shown in Fig.4.14, to study the complete behavior of this example, the simulation results are presented both for local field distributions and the global quantities as the impedance. In addition, measurements are provided to validate the simulation results from the Darwin model.

Figure 4 .

 4 Figure 4.14: Real device.

  In this sense, we are going to model the PCB framed in the black box in Fig.4.14, we do not take into account the transmission lines. The PCB model is composed of a copper loop conductor as mentioned above, shown in green in Fig.4.15, having a thickness of 35 µm and modeled in 2-D with the shell elements. The conductor is placed on a dielectric (epoxy) with height of 0.4 mm.

Figure

  Figure 4.15: PCB model.

Figure 4 .

 4 Figure 4.16: The imposed voltage sections.

Figure 4 .

 4 Figure 4.17: Considered mesh.

Figure 4 .

 4 Figure 4.19: Magnetic and electrical energies (J) vs frequency.

  pcb + L c 52.6 = L pcb + L ′ c Then, by subtracting the two equations we obtain L c -L ′ c = 75.4. (4.23) By these equations (4.23) and (4.22), we get the cable inductance L c = 111.75 nH.

Figure 4 .

 4 Figure 4.20: Magnitude of the scalar electric potential (V) in the whole domain.

  Figure 4.21: Distribution of the current density J (A.m -2 ) in the track for different frequencies.

Figure 4 .

 4 Figure 4.26: Considered model.

Figure 4 .

 4 Figure 4.27: The imposed voltage sections.

  of the impedance versus the frequencyThe evolution of the impedance as a function of the frequency obtained from the Darwin model as well as the phase are presented in Fig.4.31. The measurement results is not provided for this example. However, the capacitive-inductive coupled effects when the frequency increases can be illustrated. At low frequencies, the electromagnetic effects are decoupled. As the frequency increases, the resistiveinductive effects will be coupled, in particular, the interval [10 3 , 2 × 10 8 ] Hz. This zone corresponds to an eddy-current problem. For f > 200 MHz, the capacitive effects can be observed which indicates the importance to use the Darwin model.

Figure 4 .

 4 Figure 4.31: Modulus of impedance and the phase computed with respect to frequencies.

  In the second chapter, a state of the art on solvers linear partial differential equations is presented. Two families of linear solvers are presented in the first part: direct methods and iterative methods. In the second part, a detailed comparison considering different formulations of the Darwin model and different solvers is presented. Three Darwin formulations, in particular, M1, M2, and M3 are studied. M1 is the classic formulation of the Darwin model where no gauge technique is added.M2 is the Coulomb gauged formulation for the Darwin model. This formulation is obtained by introducing a Lagrange multiplier and rewriting the current continuity equation using the Coulomb gauge technique, which ensures the symmetry of the matrix without additional regularizations. This fact has a great effect on solving using the direct methods, where one can greatly reduce the computational time as well as the memory cost. The resulting finite element matrix derived from the M2

  The example is made for both cases. A case where the conductive track is modeled by the shell elements and the other case it is modeled by prisms. The results of comparison show a good agreement. In the last part, two examples of PCBs are modeled. A single-sided one is studied whose measurement results are provided and shown a good agreement with the results of simulations and a double-sided one is considered in order to exhibit the coupled capacitive-inductive effects and the coupling of the 2-D and 3-D conductive sub-domains. This 3-D/2-D coupling verifies the continuity of the variables (A, φ) on the interfaces of fine structures and the neighboring domains.

  electric scalar potential φ can be derived from the Maxwell equation(B.1a) such asE = -grad φ -V s grad α. (B.2) B.2. ELECTROSTATIC MODEL Then replacing the expression of E (B.2) in (B.1b), the electric scalar potential formulation for the electrokinetic model reads div (σgrad φ) = -div (σV s grad α) in Ω c , (B.3a) φ = k on Γ E . (B.3b) Let φ ′ be the test function defined on Ω. The weak formulation associated with the problem (B.3) is written Find φ ∈ H 1 (Ω) such as (σgrad φ, grad φ ′ ) Ωc = -(σV s grad α, grad φ ′ ) Ωc , ∀φ ′ ∈ H 1 Γ E (Ω). (B.4) Let φ h the approximate electric scalar potential of φ defined in the space W 0 (Ω h ) ⊂ H 1 (Ω). The weak formulation (B.4) in W 0 (Ω h ) is written Find φ h ∈ W 0 (Ω h ) such as (σgrad φ h , grad φ ′ h ) Ωc = -(σV s grad α h , grad φ ′ h ) Ωc , ∀φ ′ h ∈ W 0 (Ω h ). (B.5)

E

  electric flux density D and the electric field E. Then, the Maxwell equations (1× n | Γ E = 0, (B.7c) D • n | Γ D = 0. (B.7d) C Discrete variational formulations in space C.1 Space discretization of the EQS problem The initial form of EQS problem (1.43), described in section 1.2.2, constitues the strong formulation in φ. The weak formulation (1.86) is deduced using the Green formulas applying to the strong formulation.

  replace φ ′ h by w 0 j a vector of bases in W 0 (Ω h ), The equation (C.1) reads C.2. SPACE DISCRETIZATION OF THE MQS PROBLEM Find (φ 1 , ..., φ Nn ) ∈ R such as the problem are now carried by the scalars φ i , i = (1, ..., N n ). So, let X φ ∈ R Nn , the component vector of φ i . The problem comes down to solving the matrix systemGX φ + L Ẋφ = 0 (C.4)where G and L are the same matrices defined in (1.109) and (1.110).

  APPENDIX C. DISCRETE VARIATIONAL FORMULATIONS IN SPACEFind (A 1 , ..., A Ne ) × (φ 1 , ..., φ N n ) ∈ R × R such as ..., N n .The unknowns of the problem are carried by the scalars φ i , i = (1, ..., N n ) and the vectors A i , i = (1, ..., N e ). So, let X φ ∈ R Nn and X A ∈ R Ne , the component vector of φ i and A i respectively. And the problem comes down to solving the matrix system C 1 , G, S, and C 2 are the same matrices defined by in (1.107), (1.111),(1.109), (1.108), and (1.112).D Discrete variational formulations in timeD.1 Time discretization of the EQS problemAt time step i+1, by applying the implicit Euler scheme (1.118) to (C.4), the matrix system of the EQS model reads:

  φ) + div (σV s grad α) + div (ε ∂V s ∂t grad α) = 0, Ωc σ(grad φ + V s grad α) • grad α dΩ c + Ω ∂ ∂t ε(grad φ + V s grad α) • grad α dΩ = -I.

  8) where M, C, S, and G are the same matrices defined by in (1.107), (1.111), (1.108), and (1.109).

F

  Functional spaces and variational formulations in the shell structuresF.1 Functional spaces of differential operatorsPartial differential equations are solved in Ω s ⊂ Ω domain where Ω s is the domain containing the 2-D meshed shell surfaces and Ω is the study domain. The spaces of the differential operators are defined such that grad s , curl s , and div s belong toL 2 (Ω s ) or L 2 (Ω s ): H(grad s , Ω) = {u ∈ L 2 (Ω s ) | grad s u ∈ L 2 (Ω s )} = H 1 (Ω s ), (F.1) H(curl s , Ω) = {u ∈ L 2 (Ω s ) | curl s u ∈ L 2 (Ω s )}, (F.2)

  grad ∆ is the gradient of a triangular nodal function. The gradient is computed in the same way on the other nodes.

  1 grad s λ 2 -λ 2 grad s λ 1 )β, w 1 13 = (λ 1 grad s λ 3 -λ 3 grad s λ 1 )β, w 1 23 = (λ 2 grad s λ 3 -λ 3 grad s λ 2 )β.(G.4)

  the ζ-direction as shown in Fig. G.1.

  Figure G.1: Transformation of the real element to the reference element.

  the Jacobian matrix and e is the thickness of the thin structure. Bydenoting that J ∆ =   x 2 -x 1 x 3 -x 1 y 2 -y 1 y 3 -y 1  is the Jacobian matrix for a triangular element. The Jacobian matrix J is written as follows

  t ) -1 grad ∆ ω 0 1 β -λ 1 /e   (G.11)where |.| is the determinant of a matrix. In general, we write in the upper triangle element(J t ) -1 grad ω0 i ) -1 grad ∆ ω 0 1 β (J t ) -1 grad ∆ ω 0 2 β (J t ) -1 grad ∆ ω

  integral of the term gradgrad in the shell element is computed asΩ grad ω0 • grad ω0 dΩ = |J| Ω ((J t ) -1 grad ω 0 ) t (J t ) -1 grad ω 0 dΩ = 

1.2.2.2 Magneto-quasistatic model

					1.2. DIFFERENT MODELS
	When the displacement current represented by negligible, the MQS problem can be derived as follows ∂D in Ampere's law (1.5) becomes ∂t
	curl E = -	∂B ∂t	,	(1.44)
	curl H = J,		(1.45)
	div B = 0,		(1.46)
	div J = 0,		(1.47)
	Then replacing the expression of E (1.42) in (1.41a)-(1.41d), the electric scalar
	potential formulation for the EQS model reads
	div (σgrad φ +	∂ ∂t	(εgrad φ)) = 0 on Ω,	(1.43a)
				φ = k on Γ E	(1.43b)
	with k a constant. The scalar potential φ, defined by relation (1.42), is not unique.
	It is therefore necessary to impose a constraint on φ to guarantee the uniqueness
	of the solution. In practice, the uniqueness is obtained thanks to the boundary
	condition (1.43b), by imposing a constant value of the potential φ on Γ E .

Table 2

 2

	.1.
	1) M1 with BiCGSTAB
	2) M1 with BiCGSTAB + Split-Jacobi preconditioner
	3) M1 with BiCGSTAB + Tree Gauge + Split-Jacobi preconditioner
	4) M1 with Mumps

5) M1 with Mumps + Tree Gauge 6) M2 with Mumps 7) M3 with BiCGSTAB + Split-Jacobi preconditioner

Table 2.1: Different cases for the studied formulations.

  Ωs= -(σV s grad s α, A ′ ) Ωs -∂V s ∂t (εgrad s α), A ′

	1 µ	curl s A, curl s A ′	Ωs	+ σ	∂A ∂t	, A ′	Ωs	+ (σgrad s φ, A ′ ) Ωs
		+	∂ ∂t	εgrad s φ, A ′				Ωs	,
	σ	∂A ∂t	, grad s φ ′
									(F.6)

∈ H 1 (Ω s ) ⇒ grad s u ∈ H(curl s , Ω s ), (F.4) ∀u ∈ H(curl s , Ω s ) ⇒ curl s u ∈ H(div s , Ω s ), (F.5) ∀u ∈ H(div s , Ω s ) ⇒ div s u ∈ L 2 (Ω s ). Ωs + (σgrad s φ, grad s φ ′ ) Ωs + ∂ ∂t εgrad s φ, grad s φ ′ Ωs = -(σV s (grad s α), grad s φ ′ ) Ωs -∂V s ∂t (εgrad s α), grad s φ ′ Ωs , ∀(A ′ , φ ′ ) ∈ (H Γ B (curl s , Ω s ) × H 1 Γ E (Ω s )). (F.7)

And recalling the functional spaces taking into account the boundary conditions:

H Γ (grad s , Ω s ) = {u ∈ H(grad s , Ω s ), u = cst on Γ}, (F.8)

H Γ (curl s , Ω s ) = {u ∈ H(curl s , Ω s ), u × n = 0 on Γ}, (F.9)

H Γ (div s , Ω s ) = {u ∈ H(div s , Ω s ), u • n = 0 on Γ}. (F.10)

  F.3. DISCRETE VARIATIONAL FORMULATION Find (A h , φ h ) ∈ W 1 (Ω s h ) × W 0 (Ω s h ) such as 1 µ curl s A h , curl s A ′

					h	Ωs h	+ σ	∂A h ∂t	, A ′ h	Ωs h	+ (σgrad s φ h , A ′ h ) Ωs h
	+ εgrad s	∂φ h ∂t	, A ′ h	Ωs h	= -(σV s grad s α, A ′ h ) Ωs h -	∂V s ∂t	(εgrad s α), A ′ h	Ωs h	,
	σ	∂A h ∂t	, grad s φ ′ h	Ωs h	+ (σgrad s φ h , grad s φ ′ h ) Ωs h	+ εgrad s	∂φ h ∂t	, grad s φ ′ h	Ωs h
	= -(σ(V s grad s α), grad s φ ′ h ) Ωs h	-	∂ ∂t	(εgrad s α), grad s φ ′ h	Ωs h	,
	∀(A ′ h , φ ′

h ) ∈ W 1 (Ω s h ) × W 0 (Ω s h ). (F

.11)

  ∀i, j = 1, 2, 3, where H = (J t ) -t (J t ) -1 and |J|= |J ∆ |×e. By following the same steps, we can compute the other terms of the weak formulation obtained from (4.16) in a shell element.

	|J ∆ | |J ∆ |	e 3 ∆ e 6 ∆	(grad ∆ ω 0 i ) t Hgrad ∆ ω 0 j +|J ∆ | 1 e ∆ λ i λ j (grad ∆ ω 0 i ′ ) t Hgrad ∆ ω 0 j -|J ∆ | 1 e ∆ λ j λ i	|J ∆ | |J ∆ |	e 6 ∆ e 3 ∆	(grad ∆ ω 0 i ) t Hgrad ∆ ω 0 j ′ -|J ∆ | 1 e ∆ λ i λ j (grad ∆ ω 0 i ′ ) t Hgrad ∆ ω 0 j ′ +|J ∆ | 1 λ i λ j e ∆	              (G.14)

1.4. ASSOCIATED EXCITATION SOURCE

https://www.mouser.com/catalog/supplier/library/avxshortform.pdf

The presented Ph.D. thesis has been done in the L2EP laboratory (Laboratory of Electrical Engineering and Power electronics of Lille), in the Numerical Tools and Methods team which is co-funded by the «Hauts de France» region within the framework of FEDER (Fonds Européen de Développement Régional) and by «EDF R&Dparticipated directly or indirectly in the accomplishment of this work.

Frequency study

In the previous part, a time study in the linear and non-linear cases is made. In this part, a frequency study is presented in order to describe the behavior of the surge arrester in large frequency range. To compute the impedance of the surge arrester, it is not worth taking into account the non-linearity characteristic, a model with a linear electric constitutive law is sufficient.

In Fig. 3.13, the impedance modulus as well as the phase versus the frequency are shown. At low frequencies, in the interval [1 : 10 3 ] Hz, the surge arrester has a resistive effect, we note that the modulus of impedance is very high which means and verifies that a low current passes through the ZnO ceramic component. In the interval [10 3 : 10 7 ] Hz, a strong drop in impedance due to the presence of capacitive effects which become important. After 10 MHz, the impedance tends to be zero which means that the conductive effect of the surge arrester becomes significant.

Besides, the surge arrester has a strong capacitive-resistive coupling. However, when the excitation frequency increases, the through current becomes more and more important.

Evolution of the impedance versus the frequency

The evolution of the impedance modulus and its phase according to the frequency obtained from EQS, MQS, and Darwin models are presented in Fig. 3.30. At low frequencies, the impedance Z corresponds to the DC resistance of the winding. Then,

Time Study

In this section, we study the same example in the time domain in order to present the evolution of the current flowing out the terminal of the conductive track.

Electric potential and current density distribution

The electromagnetic fields are presented at low frequency for f = 1 Hz, the 3-D/2-D coupling between the different domains of the application is clearly shown. In

Appendices

A Ostrogradski and Green theorems

A.1 Ostrogradski's theorem

Suppose that K be a closed surface in a volume V with the boundary element dK oriented to the outside of V .

This theorem is also known under the divergence theorem which will be used for the proof of Green's formulas.

A.2 Green formulas

Let us consider Ω ⊂ R 3 be the domain of study and Γ its border. The Green formulas are defined such as

Proof. Let us consider w ∈ H 1 (Ω) and ∀v ∈ H 1 (Ω). Following the notation (1.64)

According to Ostrogradski's theorem (A.1), the term Ω div(vw) dΩ can be written:

Then, the equation (A.4) is written as follows:

We follow the same procedure for the proof of equation (A.3) and by using the fact

The electric scalar potential φ can be derived from the Maxwell equation (B.7a)

such as

Then replacing the expression of E (B.2) in (B.7b) and supposing ρ = 0, the electric scalar potential formulation for the electrostatic model reads div (εgrad φ) = -div (εV s grad α) in Ω, (B.8a)

Let φ ′ be the test function defined on Ω. The weak formulation associated with the problem (B.8) is written

Let φ h the approximate electric scalar potential of φ defined in the space W 0 (Ω h ) ⊂ H 1 (Ω). The weak formulation (B.9) in W 0 (Ω h ) is written

The problem comes down to solving the matrix system

where L is the same matrix defined in (1.110).

B.3 Magnetostatic model

In the case of a magnetostatic problem, we seek to determine the magnetic field H and the magnetic flux density B within a system subjected to an excitation. This excitation can come from the presence of permanent magnets B r or from a current flowing in an inductor modeled by an equivalent volume with a uniform distribution

Ω). (B.15)

Let A h the approximate magnetic vector potential of A defined in the space

The problem comes down to solving the matrix system

where M is the same matrix defined in (1.107) and F is a vector of dimension R Ne defined by

F.2 Continuous variational formulation

) be a test functions for the first and second equation of (4.16), respectively. Then, by applying the relations of Green on (4.16), the variational form A -φ of the Darwin model in the thin structure domain is written as: Abstract In recent years, the modeling of the magnetic and electrical components has aroused much interest in the scientific research. A magnetodynamic model is sufficient to describe the behavior of the electrical machines at low frequencies, but with the appearance of the power electronics, the machines are subjected to high frequency voltages which require the modeling of the insulators due to an accelerated aging that will be exposed. This thesis is carried out within the framework of LAMEL (joint laboratory between L2EP and EDF R&D) in order to develop a software platform namely, code carmel. This code is able to model the electrical machines in 3-D in an industrial context. This is a magnetostatic or magnetodynamic problem. Hence, the electrical fields in the dielectrics are not taken into account. The thesis is focused on quantifying and calculating the electric field in the non-conductive domain. Therefore, we are invited to implement potential formulations adapted to simultaneously calculate the electric and magnetic fields, such as the Darwin model which is able to capture the coupled capacitive-inductive effects and neglecting the wave propagation. In addition, the electrostatics and electro-quasistatic models are among the known models which are able to model the capacitive effects. Different industrial applications have been presented in this manuscript in order to validate the simulation results obtained by the Darwin model with the measurement results. Further, a Darwin thin structure problem is solved by using the shell element method, applied to both node and edge elements, in order to model the circuit boards at high frequencies. These elements are derived from the degeneration of Whitney's prismatic elements and they can be applied easily to a formulation to solve a thin shell. Keywords : Electromagnetic fields; Finite element method; Printed circuit boards; Quasistatic models; Resistive, capacitive, and inductive effects; Shell element method.

Résumé

Dans les dernières années, la modélisation des composants magnétiques et électriques suscite beaucoup d'intérêt dans la recherche scientifique. Un modèle magnétodynamique est suffisant pour décrire le comportement des machines électriques dans les basses fréquences, mais, avec l'apparition de l'électronique de puissance, les machines sont soumises à des tensions hautes fréquence, cela nécessite une modélisation des isolants en raison du vieillissement accéléré qui sera exposé. Cette thèse a été effectuée dans le cadre du LAMEL (laboratoire commun entre le L2EP et EDF R&D) avec pour but de développer une plateforme logicielle, code carmel. Ce code est capable de modéliser des machines électriques en 3-D dans un contexte industriel. Il s'agit de problème magnétostatique ou magnétodynamique. D'où, les champs électriques dans les diélectriques ne sont pas pris en compte. La thèse est focalisée à quantifier et calculer le champ électrique dans les milieux non conducteurs. Donc, nous sommes invités à mettre en oeuvre des formulations en potentiels adaptées pour calculer simultanément les champs électriques et magnétiques, telle que le modèle de Darwin qui est capable à capturer les effets capacitifs-inductifs couplés et en négligeant la propagation d'ondes. En addition, l'électrostatique et l'électro-quasistatique sont parmi les modèles connus qui sont capables à modéliser les effets capacitifs. Différentes applications industrielles ont été présenté dans ce manuscrit afin de valider les résultats de simulation obtenue par le modèle de Darwin avec les résultats de mesures. En outre, un problème de structure mince de Darwin est résolu en utilisant la méthode des éléments coques, appliquée à la fois aux éléments de noeuds et d'arêtes, afin de modéliser des circuits imprimés à des fréquences élevées. Ces éléments sont dérivés de la dégénérescence des éléments prismatiques de Whitney et elles peuvent être appliqués facilement à une formulation pour résoudre une coque mince.

Mots clés : Champs électromagnétiques ; Circuit imprimés ; Effets résistifs, capacitifs, et inductifs ; Méthode des éléments finis ; Méthode des éléments coques ; Modèles quasistatiques.