I Introduction

Français:

Le chapitre suivant introduit le sujet du chauffage par induction. Il donne un aperçu de la structure du document et des idées principales. On introduit le problème du chauffage par induction à l'égard de son application industrielle concernant un vilebrequin pour l'industrie automobile. Ensuite, on donne un bref aperçu des aspects métallurgiques et des phénomènes physiques qui ont lieu pendant le traitement thermique.

English:

The following chapter introduces the topic of induction heating. It gives an overview of the structure of the document and the main ideas of this work. It introduces the problem of industrial induction heating, with respect to its application to the surface heat treatment of an automotive crankshaft. Afterwards, a brief overview is given of the metallurgical aspects and the physical phenomena that are taking place during the heat treatment procedure.

I.1 Structure of this document

The current chapter gives an overview of the induction heating process for complex geometrical parts. The main example in this work deals with an automotive crankshaft in the context of the industrial research project OPTIPRO-INDUX. This chapter outlines the metallurgical aspects of heat treatment, as well as the major physical phenomena, which affect the heat treating procedure. The physical phenomena, as well as the metallurgical aspects are very well known. The metallurgical aspects are given in great detail in [START_REF] Barralis | Élaboration, structure-propriétés[END_REF], which has been referred to extensively in section I.4. Both the equilibrium diagram, shown in figure I.5, as well as the transformation diagram, presented in figure I.6, have been reproduced from the above-mentioned source. The experimental data for the given materials, as well as the main ideas of the theoretical concepts regarding induction heating in section I.5 can be found in [START_REF] Rudnev | Handbook of induction heating[END_REF]. These sources will not be cited again in these subsections, for reasons of brevity and readability. Additional sources complementing the former are identified.

Chapter II details the numerical model of the induction heat treatment process, including the fully transient eddy current model, the introduction of voltage potentials on closed inductor domains, as well as the heat diffusion model. It gives the reasoning behind the choice of the weak coupling procedure between electromagnetic and heat transfer computation and specifies the implementation of this coupling procedure.

The numerical model of the induction heating process features a coupled multi physics model that results in a linear system of equations, involving millions of unknowns that has to be solved in parallel with efficient numerical solvers. Chapter III presents the auxiliary space Maxwell multigrid method that is implemented for the industrial test case presented in this work. It is an algebraic multigrid preconditioner that decreases the residual error very efficiently in each iteration using suitable finite element spaces and projection operators. The theoretical background, as well as details regarding the implementation of this parallel preconditioner are given. This is followed by an application test case, including a demonstration of the advantages and shortcomings of several preconditioner combinations.

Chapter IV details a new method for the introduction of the relative movement between inductor and crankshaft, based on a discrete level set function approach. It includes the application of the method with respect to the inductor rotation around an automotive crankshaft for a large scale test case.

Chapter V gives several applications of the induction heating model for different test cases, followed by the numerical treatment of an industrial crankshaft. The crankshaft model is provided by industrial partner PSA Peugeot Citroën and includes an inductor provided by EFD Induction. The material for the crankshaft is provided by ASCOMETAL.

Finally, chapter VI sums up the key results of this work and gives an outlook on directions for possible future research.

I.2 Induction heat treatment

The group of industrially produced metallic workpieces is very diverse. Some workpieces used in the automotive sector include springs, wires, camshafts, brake disks, screws or crankshafts. Their production necessitates several distinct steps. In most cases, the initial raw part is formed, brought into shape and in many cases made to connect with other raw parts to form a new part. Afterwards or sometimes in between, the product undergoes further treatment to change its appearance or material properties. These distinct varieties of the production cycle are achieved by utilising a wide variety of manufacturing methods. The initial shaping can be realised by rolling or forging, whereas additional connections can be made by glueing, cold and hot welding or soldering. The appearance of the workpiece is mostly changed using mechanical techniques like filing, grinding, honing or shaving or chemical treatments like galvanisation.

The material behaviour, like the hardness or the ductility is influenced by the material composition, as well as by the microstructure evolution of the material. The internal distribution of the material composition, as well as the related microstructure depend on the chosen production process, the amount of applied deformation energy, the duration of each procedure, as well as the order of each subsequent production step. The application requirements of the finished product are as diverse as the application areas. Workpieces that are stretched, like springs should remain ductile during their planned use, yet they must retain their initial shape, while unloaded. Brake disks are used to transform kinetic energy into heat using friction. These workpieces must be hard, in order to prevent surface degradation and to prevent fatigue and crack growth, while ensuring a consistent connecting area and good heat conducting properties in order to extract the generated heat from the application area. Screws are used as fasteners to connect different workpieces. They are under a constant tensile stress. It is therefore very important that they resist the effects of creep. Yet, screws must also be able to withstand shocks under changing load cycles with imposed tensile and shear stresses. Crankshafts or camshafts are used to translate rotational motion into vertical motions or vice versa. They are supported by bearings. These components must endure changing loading cycles, shocks and temperature differences. Their cores must show a ductile behaviour. In contrast, the surface connection between the shafts and the supporting bearings or different components like valves for camshafts necessitate a hardened surface, such that the effects of wear are reduced. This very diverse set of requirements is often achieved by utilising a heat treatment procedure.

The mechanical behaviour of steel is related to its internal microstructure. An introduction to the general properties of metals and microstructure evolution can be found in [START_REF] Barralis | Élaboration, structure-propriétés[END_REF]. The ductility and general softness of steel is related to the grain size and shape of the microstructure. Ductility can be increased by annealing. For this procedure the temperature of the material is increased and held above a point at which diffusion happens more easily. First, internal dislocations are reduced and the internal grain energy is minimised in the recovery stage. Afterwards, recrystallisation and grain growth lead to a homogenisation and an increase in grain size. The large grains that are bounded by smooth borders result in a softer material. Hardening by introducing a martensitic phase is another heat treating procedure. It is achieved by elevating the temperature of the material up to a point, where the material configuration changes into the austenitic phase. Afterwards, the material is rapidly cooled by quenching to force a displacive, i.e. diffusionless transformation into its martensitic form. The large quantity of dislocations results in a high material hardness.

The heat energy, which is necessary for elevating the process temperature can be introduced using a wide array of methods, like convective heating using furnaces, direct electrical resistive heating or induction heating. Convective heating utilises heat diffusion from an outside heating source to change the temperature distribution inside the workpiece. This results in a distribution of heating energy inside of the workpiece that is mostly uniform, with a temperature gradient starting from the surface, extending up to the core of the workpiece. When the process time is increased the temperature gradient is minimised, such that the workpiece has a homogeneous temperature. Joule heating is a direct form of heating that is based on prescribing currents by generating a voltage potential on opposing contact surfaces of the conductive workpiece domain. The internal current generates heat, due to Ohm's law. The heating efficiency is very high, but the process itself is difficult to apply to workpieces of complex shape, because the current can not be applied locally. It can only be applied in a global manner, such that it is difficult to focus the current only on certain features of the geometry, like the surface. Electromagnetic induction heating follows the same approach as resistive Joule heating in that it induces currents into the workpiece, which then in turn generate heat. The distinction of electromagnetic induction heating is that the voltage potential is not imposed using direct contact. Rather, an alternating current is imposed on an inductor that is placed near the workpiece surface. This changing electric current creates a changing magnetic field, which induces eddy currents inside the conducting domain of the workpiece. The skin effect will lead to a concentration of eddy currents near the surface of the conducting domain. The amount and penetration depth of the applied heating energy can be influenced by the shape and magnitude of the applied loading. This permits the contactless heating of conductive workpieces of complex geometrical shape. The process parameters can be quickly adapted according to the needs of the producer. The high heating rates result in a cost effective procedure. One of the major advantages of induction heating is that alternating currents tend to concentrate near the surface of conducting domains, due to the skin effect. This allows to concentrate the heating energy in the regions of interest, without affecting the material core. It is this feature that makes induction heating a very suitable method for use in surface hardening applications.

The goal of this work is to introduce a numerical model for the surface hardening procedure by electromagnetic induction for an automotive crankshaft. These workpieces are commonly formed by forging, followed by a preheating step to elevate the workpiece to a uniform initial temperature. Afterwards, heating sources are applied to the surfaces that will later support bearings, using induced eddy currents. After a material and process dependent time the workpiece is rapidly cooled down in a water or oil bath to form the martensite in the region of the contact surface. A possible tempering step can be added, in order to reduce the large internal stresses that are generated by the hardening procedure. The steps for this procedure are comparable to annealing. After quenching the workpiece is heated to an elevated temperature, which results in a reduction of hardness and internal stress.

I.3 Industrial induction heating of a crankshaft

Crankshafts are heat treated in the regions that will support bearings. In the beginning, the steel workpiece is clamped to a frame, as can be seen in figures I.2 and I.3. This allows to rotate the crankshaft, with respect to the inductors, so that the heating can be generated as evenly as possible. Figure I.2 shows the inductor assembly. It is housed in a casing that features cooling pipes and sometimes measuring devices, which can register the temperature. The inductors, like the crankshaft, will heat up during the heat treatment. Therefore, the inductors are liquid cooled during the whole procedure. The heated workpiece is finally quenched after a certain time. This can be achieved, either by placing the workpiece in a water bath, as can be seen in figure I.3 or by spray cooling, as is depicted for a heat treated workpiece in figure I.4.

I.4 Metallurgical aspects of hardening

Hardening of steel is accomplished by changing the microstructure configuration of the material. The workpiece is usually pre-heated up to a certain temperature T p . Afterwards heat is distributed to the regions of interest to further increase the temperature T up to a temperature T γ , which is the point where an austenitic phase γ forms. In general, T γ is chosen, such that it is well above the transition zone of each distinct phase to ensure that the regions of interest are well in austenitic configuration. Yet, it must not be too high in the heat affected zone, so that the workpiece does not remain in the austenitic region for an extended period of time during the cooling phase. The precondition for the martensitic phase transformation is the stable existence of ferrite α and carbide Fe 3 C at ambient temperature T A . As an example, figure I.5 shows a simplified equilibrium diagram for a non-alloyed steel in the low-carbon content range. The practical temperature range for T γ for non-alloyed hypoeutectoid and hypereutectoid steels is indicated in red and green. For hypereutectoid steel it is visible that the martensitic phase transformation can either start from a phase configuration of pure γ or a mix of γ and carbides. This temperature range is shown in blue.

The martensitic phase transformation results from a displacive transformation of austenite γ into martensite M . This change is rapid, nearly instantaneous and diffusionless. The dissolved carbon in γ can not exit the crystal structure, which results in the new martensitic configuration M , including the carbon in a desequilibrium condition. The crystallographic microstructure of martensite is built from crystals including carbon, preferably on the crystal borders. This results in a large number of dislocations. The material is in a constant stress/strain Table I.1: Composition of hypoeutectoid steel 41Cr4, as given by [START_REF] Bohemen | The kinetics of bainite and martensite formation in steels during cooling[END_REF] configuration, which substantially increases the material hardness. The quenching process is highly dependent on the cooling rate of the workpiece. Figure I.6 shows a simplified form of the continuous cooling transformation (CCT) diagram for an alloyed hypoeutectoid steel of type 41Cr4, with the composition given in table I.1. The diagram shows the microstructure evolution, with respect to the cooling time for material samples that have been austenitised for 30 min at 850 • C. The region, where austenite transforms into martensite, is shown in light-blue. Above T Ms exists the unstable γ-phase that becomes stable at a certain temperature T ac1 . The stable area is shaded in light red. The light-green shaded area identifies the region, where ferrite forms, whereas the dark-green and darkblue shaded areas indicate the regions, where the phase transforms into a perlitic and a bainitic structure. The curves ①,② and ③ show the microstructure transformations for samples that have been cooled at three different cooling rates. Curve ① is rapidly cooled, such that the fraction of martensite reaches y M = 1, after the quenching is finished. Curve ② is cooled more slowly, such that it traverses both the perlitic, as well as the bainitic regions, before entering the region, where austenite transforms into martensite. At ambient temperature, the resulting material composition will include a low fraction of martensite y M ≪ 1. The remaining fractions will consist of a ferritic, a perlitic, as well as a bainitic microstructure. Curve ③ shows the transformation of a very slowly cooled sample. It never reaches the region, where austenite transforms into martensite, such that the final fraction of martensite at ambient temperature will be y M = 0. The final microstructure will be a composition of a ferritic, as well as a perlitic phase.

The martensitic phase transformation starts at a given temperature T Ms , such that the evolution of martensitic fraction y M and austenitic fraction y γ can be written with respect to a lower temperature T l using an empirical law by [START_REF] Koistinen | A general equation prescribing the extent of the austenite-martensite transformation in pure iron-carbon alloys and plain carbon steels[END_REF], as y M = 1exp(-k(T Ms -T l ) n ).

(I.1)

The relative amount of martensite can be given for the temperature change

∆T M sl = T Ms -T l , (I.2)
with respect to the fraction of remaining austenite y γ as

∆T M sl = n  ln(y γ ) -k . (I.
3)

The starting temperature for the martensitic transformation T Ms , as well as the final temperature of the martensitic transformation strongly depend on the alloys of the steel. Physical elements like Cr, Mo, Si, Ni or Mn can be added to the steel, in order to reduce the cooling velocity. The hardenability of a material is related the cooling time, for which the material goes directly from the austenitic phase to a mix of austenitic and martensitic phases. E.g. the steel 36NiCrMo16 forms a final microstructure decomposition of pure martensite for a sample with a diameter of 10 mm, when allowed to be cooled by the surrounding air. In contrast, for a sample of a steel of type 2C60 the same final microstructure can only be reached by quenching in a water bath. The former steel possesses, therefore, a greater hardenability as 2C60. E.g. [START_REF] Bohemen | The kinetics of bainite and martensite formation in steels during cooling[END_REF] states that the temperature T Ms for the hypoeutectoid steel 41Cr4 lies approximately at T Ms = 304 • C and that both empirical parameters of equation (I.4) can be given as k = 0.016 and n = 1, such that 50 % of martensite will have formed for a reduction of temperature of ∆T M sl = ln(0.5) -0.016 ≈ 43.32 • C, (I.4) whereas 90 % of martensite will have formed after a reduction of temperature of

∆T M sl = ln(0.1) -0.016 ≈ 143.91 • C. (I.5)
The final martensitic steel is very hard and can be very brittle. A tempering step can be added to the heat treatment procedure to reduce the large internal stresses in the workpiece. Figure I.7 shows two hypothetical heating cycles for a production process. First, the workpiece is pre-heated, afterwards heat is applied to reach T γ . After a process dependent time t γ the workpiece is quenched to force the martensitic transformation. The first process, shown in red, reheats the workpiece to an elevated tempering temperature T t 1 , before the ambient temperature T A is reached. In the second process, which is shown in blue, the workpiece is fully cooled to T A before the tempering step begins. A second tempering step with a lower temperature T t 2 is added to further change the material configuration. Alternating currents tend to concentrate near the surface of the conductor, due to the changing magnetic field. This magnetic field generates a changing electric field, which leads to a concentration of eddy currents towards the surface of the domain.

It results in a strong concentration of current near the surface, which is far greater than at the centre of the workpiece, for which the current density vanishes, due to current cancellation effects. An increase of the frequency for the applied current leads to a higher concentration of eddy currents near the surface. The current density distribution is usually approximated using a decaying function, such that the current density per unit area at the surface I s decays exponentially as

I = I s exp(- x δ
), (I.6) such that the current density per unit area I inside of the domain decays up to

I = I s e (I.7)
for a distance x = δ. The value δ is the penetration depth of the material. The assumption is that the conductive domain has an infinite dimension, such that the surface can be set to x = 0, whereas the centre lies in infinite distance. The relative amount of current that is distributed inside the workpiece, between x = 0 and x = δ, can then be approximated as

 δ 0 Idx  ∞ 0 Idx ≃  δ 0 exp(-x)dx  ∞ 0 exp(-x)dx
≈ 63.21 %.

(I.8)

The heating power density Q em is related to the current density J and the electrical resistivity ρ as

Q em = ρJ 2 . (I.9)
It follows that the relative amount of heating power that is distributed inside the workpiece, between the surface and x = δ can, therefore, be approximated as

 δ 0 Q em dx  ∞ 0 Q em dx ≃  δ 0 exp(-2x)dx  ∞ 0 exp(-2x)dx
≈ 86.47 %.

(I.10) Accordingly, less than 15 % of the heating power during an induction treatment is distributed further inside the workpiece domain than the penetration depth δ.

The penetration depth can be approximated for a homogeneous solid with constant resistivity and magnetic permeability µ, by µ = µ 0 µ r (I.11) using the frequency of the induced current f , such that

δ =  2ρ 2πf µ =  ρ πf µ 0 µ r =  ρ πf 4π10 -7 H m µ r ≈ 503  ρ f µ r H -1 2 m 1 2 .
(I.12)

Here, µ 0 is the magnetic constant, often called magnetic permeability in free space and µ r is the relative magnetic permeability of the conducting material. The magnetic permeability in free space µ 0 is arbitrarily defined [START_REF] Mohr | CODATA Recommended Values of the Fundamental Physical Constants[END_REF], such that two current-carrying wires of length l = 1 m, seperated by a distance d = 1 m, carrying each a current of J w = 1 A attract each other by a resulting force of F = 2 • 10 -7 N, such that

µ 0 = 4π • 10 -7 H m .
(I.13) Figure I.8 shows the skin effect for the current density per unit area for aluminium at 20 • C for a current,which is alternating with a frequency of 10 kHz. The shadowed area marks the penetration depth δ = 0.826513 mm. At this temperature it can be assumed that the electrical resistivity is ρ = 0.027 • 10 -6 Ωm, whereas the magnetic permeability is very close to the magnetic permeability of free space, such that µ r = 1. In general, the penetration depth increases with rising resistivity and progressively smaller values of the magnetic permeability. It is visible that the penetration depth is inversely proportional to the frequency of the alternating current. It can be assumed that metallic materials possess a strong temperature dependency with respect to the resistivity and relative magnetic permeability. E.g. the resistivity of stainless steel and aluminium increases with rising temperature. The relative magnetic permeability of aluminium is close to that of free space, whereas the relative magnetic permeability of stainless steel is higher, but decreases with rising temperature. At a certain material dependent temperature T C , which is called the Curie temperature, the material properties of steel change abruptly, resulting in a relative magnetic permeability of µ r = 1. Figure I.9 shows the effects of increasing temperature to the penetration depth δ for general steel workpieces. Initially, the penetration depth increases steadily with rising temperature, due to increasing resistivity and decreasing relative magnetic permeability. After reaching the Curie temperature T C the relative magnetic permeability becomes suddenly µ r = 1, which results in a large increase of the penetration depth. Afterwards, the electrical resistivity rises further with increasing temperature.

Figure I.10 shows the effects of changing frequency on the induced currents for aluminium at 20 • C for f = {1, 10, 100} kHz, which results in penetration depths of δ = {2.61366, 0.8263, 0.26137} mm (see also table 2.1-1 in [START_REF] Jürgens | Induction heating and melting[END_REF]). The current inside the domain approaches I s for lim f → 0. From equation (I.12) follows that under the assumption of equal resistivity, the penetration depth of a material with relative magnetic permeability of µ r = 100 under the influence of an alternating current with a frequency of f = 100 Hz will have the same penetration depth as a material with a relative magnetic permeability µ r = 1 under the influence of an alternating current with a frequency of 10 kHz. A large relative magnetic permeability has, therefore, the effect of pushing the induced eddy currents away from the relative centre of the material. The frequency input, as well as the duration of the heat treatment, are the defining properties for a surface hardening process. The frequency has a direct impact on the depth of the induced eddy currents, whereas the duration of the treatment influences heat diffusion in the workpiece. This observation leads to the conclusion that a homogeneous heating is accomplished using low frequencies, whereas surface heat treatment is done in the high frequency domain. In practice, it must be assumed that the workpiece material is not homogeneous, which is the case for most of the industrially used materials and that the workpiece is not of simple shape, such that equation (I.12) can only be used to approximate the penetration depth. It follows that the induced electrical current will in general not follow the path of exponential decay, due to non-linear material behaviour. This motivates the implementation of a numerical model for the coupled induction heating process.

I.5.2 Magnetic permeability and Curie temperature

The magnetic permeability µ is the constitutive relationship between the magnetic flux density B and a magnetising magnetic field H such that B = µH. (I.14) This simple relationship states that materials with high relative permeability µ r allow the easy formation of magnetic fields. All materials can be classified into three sub-categories. Dielectric materials possess a relative magnetic permeability that is smaller than 1, which means that it is more difficult to form a magnetic field in these materials than in free space. Diamagnetism is a weak effect present in all materials. It is the tendency of the material to create an opposing magnetic field to an external field. Paramagnetic materials possess a relative magnetic permeability µ r ≥ 1, whereas ferromagnetic materials possess a relative magnetic permeability µ r ≫ 1. The distinction between the latter materials is that the crystallographic structure of ferromagnetic materials is made up of Weiss domains [START_REF] Kurz | Traité des Matériaux[END_REF], as can be seen in figure I.11b. Figure I.11: Difference between paramagnetic and ferromagnetic materials These domains include regions with equal parallel magnetic moments, produced by the alignment of the magnetic dipoles inside the material. These domains can easily adjust in the direction of an applied magnetic field, whereas paramagnetic materials possess a chaotic distribution of magnetic moments, as can be seen in figure I.11a. For these materials, an alignment is possible, once the material is exposed to a magnetic field, but it has to be larger than for the ferromagnetic case.

Ferromagnetic materials show an inherent magnetisation M , such that relation (I.14) should be rewritten as

B = µ 0 (H + M ) = µ 0 (1 + M H )H = µ 0 (1 + χ)H, (I.15)
with the magnetic suscebitibility χ. Ferromagnetic materials retain a magnetic field M , which forms spontaneously [START_REF] Kurz | Traité des Matériaux[END_REF]. For these materials, the relationship between both magnetic flux density B and applied magnetic field H is highly non-linear as can be seen in figure I.12. Initially, B increases proportionally Figure I.12: Saturation effects for ferromagnetic materials to H up to a certain point M s , which is the saturation point for the specific material. It is the maximum attainable magnetic field M inside the material. The application of an outside magnetising field H leads to a reorientation and resizing of the existing Weiss domains in a ferromagnetic material. This effect, leads to a change in remaining magnetisation M r that remains, even when the outside excitation is turned off. This magnetic hysteresis effect is visualised in figure I.13. The initial starting point with zero remaining magnetisation is marked as point ①. The initial excitation using an outside magnetising field H up to saturation point ② leads to a build up of M r . Point ③ shows the reminiscence that exists without outside excitation. The remaining magnetisation vanishes, if a magnetising field is created that points in the opposite direction. The field strength that is necessary, such that M r vanishes is given as the distance between point ④ and point ① and is defined as the coercive field. A changing magnetic field will create a hysteresis loop, that goes up to the bottom saturation point ②' and follows through ③' to ④'. A rapidly changing magnetic field, like in induction heating applications, will lead to a repetitive resizing and reorientation of the existing Weiss domains inside the material. This will lead to hysteresis losses, which result in the formation of thermal energy. In standard working steels it can be assumed that this effect is negligible, in comparison to the heat generated by eddy currents inside the work piece.

Increasing material temperature leads to increased thermal motion of the atoms of the material. This has the effect that the magnetic monopoles can not align any more into Weiss domains. This means that the initial ordering of the magnetic monopoles becomes chaotic, like for the paramagnetic case. Figure I.14 shows the diminishing internal magnetisability M for rising temperature, with maximum value M 0 , in green. T C is the curie temperature, at which the material becomes paramagnetic. This change in phase occurs simultaneously with a rapid change of specific heat capacity C, which is shown in red. Let F be the Helmholtz free energy, G the Gibbs free energy and let S be the entropy of the system. [START_REF] Gallagher | Handbook of thermal analysis and calorimetry[END_REF] shows that the entropy is related to the first derivative of F for whereas the specific heat capacity is related to the second derivative

C = - ∂ 2 F (∂T ) 2 . (I.17)
The classification of Ehrenfest [START_REF] Ehrenfest | Phasenumwandlungen im ueblichen und erweiterten Sinn, classifiziert nach den entsprechenden Singularitaeten des thermodynamischen Potentiales[END_REF] states that a first order phase change is defined through the emergence of a discontinuity in at least one state variable, like the entropy (I.16). A second order phase change happens, if the discontinuity first occurs in a second derivative, like the specific heat capacity (I.17). [START_REF] Mattis | Statistical mechanics made simple[END_REF] states that the magnetisation is a state variable for magnetic materials, such that it can be defined by the first derivative of the Helmholtz free energy, for constant T (cf. [START_REF] Callaway | Quantum theory of the solid state: Part A[END_REF]), as

M = - ∂F ∂H , (I.18)
such that the magnetic susceptibility is proportional to the second derivative. It can be written with respect to B as It can be seen that C is not defined for point T pc . Standard engineering materials, like steel are multi-component systems. In these cases, the first order phase change does not happen as an isothermal process, but in a defined temperature range. Figure I.15b shows a second order phase change at the Curie temperature T C . The specific heat capacity changes rapidly in a discontinuous manner, but remains defined [START_REF] Callaway | Quantum theory of the solid state: Part A[END_REF]. The magnetisation (I.18) will continuously increase from zero below the Curie temperature to its maximum M 0 . The discontinuity exists for the magnetic susceptibility in (I.19) and the specific heat capacity C in equation I.17, which will result in a jump for both quantities.

χ = ∂ 2 G (∂B) 2 . (I.19)
For each specific heat capacity, shown in figures I.14 and I.15, the heat capacity increases from an initial starting value C 11 to C 12 at the critical phase transition point. In both of the presented cases, the new starting value C 21 for a temperature T > T C differs from C 11 . With rising temperature each specific heat capacity increases to its final value C 22 for the given temperature range. It must be noted that steel in γ-phase does not possess a Curie point. The austenitic phase is, therefore, fully paramagnetic, such that µ r = 1. The above-mentioned hysteresis effects can, therefore, not exist in austenitic steels.

II Numerical model for induction heating

Français:

English:

The following chapter introduces the physical model of electromagnetics based on Maxwell's equations, as well as the general heat diffusion model. This is followed by a derivation of a hyperbolic and a parabolic version of the electromagnetic equations. Using the classification principles for partial differential equations, in connection with the method of characteristics, it is shown that the parabolic magnetic vector potential form is more adequate to describe the induction heating problem. This section is followed by a brief introduction to the Galerkin finite element method, in connection with the curl-conforming Nédélec vector finite elements and the weak formulations of the magnetic vector potential equation, the heating equation and the voltage potential problem. This overview over discretisation methods is finalised by an introduction of stable time discretisation methods.

Finally, an algorithm is presented to impose conforming source currents on arbitrarily complex closed inductor geometries, which has been identified as prerequisite for the well-posedness of the underlying electromagnetic problem.

II.1 Physical model

II.1.1 Maxwell's equations

The fundamental equations for describing the electromagnetic phenomena are given by Maxwell's equations in integral form (see, e.g. [START_REF] Jin | The finite element method in electromagnetics[END_REF]). Maxwell's equations are four fundamental equations that are Faraday's law, the Maxwell-Ampère law, as well as the flux theorem for the electric flux density of Gauss and the condition of divergence free magnetic fields. Let dl, ds and dv be the smallest differential forms of a line, a surface area and a volume segment. Faraday's law states that the negative change of the time derivative of the magnetic flux density B of a given surface area equals the electric field E over the contour line of this surface area as

 line Edl = - d dt  surface
Bds.

(II.1)

The Maxwell-Ampère law relates the magnetising field H, expressed through a line integral, with the sum of the time derivative of the electric flux density D and the electric current density J integrated over a given area, such that (II.

3)

The former equations are taken together with Gauss's magnetic law, which states that every magnetic field must be divergence free, such that its surface integral over a closed surface must vanish as

 surface Bds = 0. (II.4)
Let x, y, z be the x-,y-and z-direction of the 3-dimensional space R 3 and e x ,e y and e z the unit vectors defining the respective coordinate system. Let ζ v be an arbitrary differentiable vector field in a given vector space V in R 3 and let ζ s be an arbitrary quantity in the scalar field S. Further, let ∂ζs ∂x , ∂ζs ∂y and ∂ζs ∂z be the spatial partial derivatives of each scalar quantity of ζ v with respect to each coordinate direction.

The operator ∇ can be defined as the gradient operator mapping from a scalar field to a vector space ∇• can be defined as the divergence operator that maps from a vector field into a scalar space

∇ : s ∈ S → v ∈ V, (II.
∇• : v ∈ V → s ∈ S, (II.7) such that ∇ • ζ = ∂ζ x ∂x + ∂ζ y ∂y + ∂ζ z ∂z . (II.8)
∇× can be defined as the curl-operator, which is an operator mapping from a vector to another vector

∇× : v ∈ V → v ∈ V, (II.9) such that ∇ × ζ = ( ∂ζ z ∂y - ∂ζ y ∂z )e x + ( ∂ζ x ∂z - ∂ζ z ∂x )e y + ( ∂ζ y ∂x - ∂ζ x ∂y )e z .
(II.10)

Stokes's theorem (see, e.g. [START_REF] Gurtin | An Introduction to Continuum Mechanics[END_REF]) states that the closed line integral of a vector field ζ v can be represented by the surface integral over any open surface, bounded by line l as

 surface ∇ × ζ v ds =  line ζ v dl. (II.11)
Gauss's integral theorem (see, e.g. [START_REF] Gurtin | An Introduction to Continuum Mechanics[END_REF]) states that the flow of a vector field ζ v over a closed surface boundary equals the volume integral of the divergence, such that

 volume ∇ • ζ v dv =  surface ζ v ds.
(II.12) Equations (II.1) and (II.2) can be rewritten in differential form (cf. the original reformulation by [START_REF] Heaviside | Electrical papers[END_REF], Art. 30, sec. 18) using Stoke's theorem (II.11) and the localisation theorem of continuum mechanics (see, e.g. [START_REF] Gurtin | An Introduction to Continuum Mechanics[END_REF]), as

∇ × E = - ∂B ∂t (II.13)
and

II.1 -PHYSICAL MODEL II -NUMERICAL MODEL ∇ × H = ∂D ∂t + J. (II.14)
Under the same assumptions, equations (II.3) and (II.4) can be rewritten using Gauss's theorem (II.12) as

∇ • D = γ. (II.15)
and

∇ • B = 0. (II.16)
The divergence of a rotational field vanishes (see, e.g. [START_REF] Girault | Finite Element Approximation of the Navier-Stokes Equations[END_REF]), such that the application of the divergence operator (II.7) to equation (II.14), with subsequent application of equation (II.15), leads to the equation of continuity for the electric current density J in relation to the electric charge density γ as

∇ • J = - ∂γ ∂t , (II.17)
which can also be written in integral form using (II.12) as

 surface Jds = - d dt  volume γdv.
(II.18) Equation (II.16) follows in the same manner from equation (II.13). The differential form needs to be completed with the constitutive relationships, in order to include the electromagnetic properties of the materials, using the electrical permittivity ϵ, the magnetic permeability µ, as well as the electrical conductivity σ, as

D =ϵE (II.19) B =µH (II.20) J =σE. (II.21) II.1.1.1 E-H-J formulation
The E-H-J-formulation can be derived after inserting the constitutive equation (II.19) into (II.14) and (II.20) into (II.13) and reordering, such that

∂ϵE ∂t -∇ × H = J ∂µH ∂t + ∇ × E = 0. (II.22) (II.23)
Subsection II.2.1.3 demonstrates that the system of equations (II.22) and (II.23) is of hyperbolic type. Therefore, it is difficult to handle numerically, due to the wave-like nature of the propagation. This is especially problematic in a global finite element setting, like the one presented in this work, due to scattering effects and the use of distinct media. The E-H-J formulation has been effectively used in the context of electromagnetic scattering applications, using explicit formulations [START_REF] Hesthaven | Nodal High-Order Methods on Unstructured Grids[END_REF], based on a decoupled explicit nodal discontinuous Galerkin scheme [START_REF] Hesthaven | Nodal discontinuous Galerkin methods: algorithms, analysis, and applications[END_REF]. In these formulations a wave transport problem is solved, including scattering and refraction effects. For the industrial setting of induction heating it can be assumed that the electromagnetic problem is in a sufficiently stationary condition, so that it possesses a quasi-stationary behaviour. Therefore, the following diffusion-like magnetic vector potential formulation is more adequate to describe the electromagnetic phenomena.

II.1.1.2 A-V formulation

The current density in (II.14) can be split into an induced current J d that depends on the changing magnetic field and an imposed source current J s , such that

J = J d + J s .
(II.24) [START_REF] Heaviside | Electrical papers[END_REF] shows that the magnetic flux density can be expressed by a magnetic vector potential A, such that

B = ∇ × A. (II.25)
After inserting (II.25) and (II.21) into Faraday's law (II.13) and utilising the linearity of the curl operator (II.9) the dependent current J d can be expressed by

J d = -σ ∂A ∂t . (II.26)
According to [START_REF] Jin | The finite element method in electromagnetics[END_REF], the source currents can either be expressed directly or by introducing a voltage potential Φ, such that E = -∇Φ.

(II.27)

Inserting equation (II.27) into the charge conservation equation (II.17) and using the constitutive equation (II.21) leads to

∇ • σ∇Φ = ∂γ ∂t .
(II.28)

The displacement currents ∂D ∂t can be neglected in induction heating (see, e.g. [START_REF] Rudnev | Handbook of induction heating[END_REF]). In this context, the Maxwell-Ampère law (II.14) can be rewritten in the A-V-formulation using (II.25) and (II.28), assuming ∂γ ∂t = 0 in the imposed source term, as

σ ∂A ∂t + ∇ × 1 µ ∇ × A + σ∇Φ = 0 ∇ • σ∇Φ = 0.
(II.29) (II.30) Subsection II.2.1.2 demonstrates that equation (II.29) is a parabolic partial differential equation, whereas (II.30) is of elliptic type. Both equations can be weakly coupled, under the assumption of negligible influence of the electromagnetic phenomena on the voltage potential of the source current, so that each equation can be treated separately in an efficient manner.

The A-V formulation is widely used in the literature, e.g [START_REF] Ren | Influence of the R.H.S. on the Convergence Behaviour of the Curl-Curl Equation[END_REF] and [START_REF] Biro | On the Use of the Magnetic Vector Potential in the Nodal and Edge Finite Element Analysis of 3D Magnetostatic Problems[END_REF] use this formulation in a global finite element setting in stationary form. [START_REF] Hiptmair | Coupled boundary-element scheme for eddy-current computation[END_REF] uses a time harmonic formulation in connection with a boundary element formulation. [START_REF] Houston | Mixed Discontinuous Galerkin Approximation of the Maxwell Operator: Non-Stabilized Formulation[END_REF] introduces a discontinuous Galerkin approach for the stationary vector potential form and [START_REF] Kolev | Parallel Auxiliary Sace AMG for H (curl) Problems. Technical report[END_REF] uses the fully transient form for a highly parallelised application involving up to approximately 80 million degrees of freedom on 1024 processor nodes.

An alternative potential formulation can be derived based on a magnetic scalar potential in connection with an electric current vector potential [START_REF] Biro | Computation of 3-D Magnetostatic Fields Using a Reduced Scalar Potential[END_REF] (see also [START_REF] Preis | Different Finite Element Formulations of 3D Magnetostatic Fields[END_REF] and [Biro et al., 1993a]). Applications to an induction hardening application are presented in [START_REF] Candeo | Multiphysics Modeling of Induction Hardening of Ring Gears for the Aerospace Industry[END_REF]. It states that with respect to a comparable A-V potential formulation the advantage of this formulation is a reduction in storage costs and solution time for the electromagnetic problem, even though the accuracy is reduced.

II.1.2 Heat transfer equations

The first law of thermodynamics states that the total energy of a system can only change due to incoming and outgoing heat fluxes, work done by or on the system or as a result of an internal change in energy due to physical phenomena (see, e.g. [START_REF] Wriggers | Nonlinear finite element methods[END_REF]). These physical phenomena include radioactive decay, heat generation due to induced eddy currents Q em or heat generation or absorption due to phase changes (see the definition of T C and the magnetic hysteresis effects in section (I.5.2)). The internal and external work, as well as the changes in kinetic energy, can be neglected in the context of induction heating. Therefore, the internal energy can be stated as a function of internal heat. The internal effects of radiation and phase changes are much weaker than the heating energy that is generated by induced eddy currents. Thus, the integral form of the energy conservation law can be given by d dt

 volume ρCT dv = -  surface q F ds +  volume Q em dv,
(II.31) using the density ρ, the specific heat capacity C, defined in equation (I.17), the induced heating energy Q em , defined in equation (I.9), as well as the surface heat flux over the enclosing domain q F , which points in the opposite direction of the surface normal. The differential form can be written, using Gauss's theorem (II.12), such that

∂ ∂t ρCT + ∇ • q F = Q em . (II.32)
Fourier's law (see, e.g. [START_REF] Mortimer | Physical Chemistry[END_REF]) relates the heat flux to the gradient of the temperature field T , using the thermal conductivity k, such that q F = -k∇T.

(II.33) Equation (II.32) can, therefore, be rewritten in the form

∂ ∂t ρCT -∇ • k∇T = Q em .
(II.34)

The heating flux q F may describe the heat loss due to temperature differences between the workpiece and the surrounding medium q Fm or the heating loss following from the emission of electromagnetic radiation q Fr as demonstrated in [START_REF] Bay | A numerical model for induction heating processes coupling electromagnetism and thermomechanics[END_REF]. The convection flux can be described using the ambient temperature T A and a convection coefficient h as

q Fm = h(T -T A ), (II.35)
whereas the radiation losses q Fr can be defined as

q Fr = εσ s (T 4 -T 4 A ), (II.36)
using the material emissivity ε and the Stefan-Boltzmann constant σ s , as defined in [START_REF] Mohr | CODATA Recommended Values of the Fundamental Physical Constants[END_REF].

In the following sections it is shown that the parabolic equation (II.34) is easier to handle numerically than equation (II.32), which is hyperbolic. The boundary conditions (II.35) and (II.36) are not introduced directly in equation (II.34). Both conditions are introduced after the introduction of the weighted residual formulation and integration by parts with application of Green's theorem, as shown in subsection II.2.2.2. [START_REF] Zwillinger | Handbook of differential equations[END_REF] states that before any given differential equation is approximated, it should be established, whether it is inherently well-posed. The procedure follows the criteria of Hadamar (see, e.g. [START_REF] Lanczos | Linear Differential Operators[END_REF]), which define that a differential equation is well-posed if:

II.2 Numerical model

II.2.1 Classification of PDEs

• The solution exists.

• The solution is unique.

• The solution is stable, in the sense that small perturbations of the input data or boundary conditions do not lead to highly different outcomes for the solution.

The first two criteria, which are existence and uniqueness depend not only on the partial differential equation itself, but on the chosen solution method. Subsection II.2.4 establishes these conditions for the chosen methods of this work. The third criterion, which is stability, also depends on the chosen solution method, but in addition it is also inherently related to the characteristic behaviour of the partial differential equation, i.e. the underlying problem that is described.

Partial differential equations (PDEs) can be divided into three subclasses showing a different numerical behaviour. The subclasses are elliptic, parabolic and hyperbolic partial differential equations. In the following a second order PDE with two independent variables will be used to describe the classification theory. The classification of higher dimensional problems follows analogously. [START_REF] Zwillinger | Handbook of differential equations[END_REF] states that any second order PDE of the above-mentioned form can be written as

A(x, t) ∂ 2 u ∂x 2 + B(x, t) ∂ 2 u ∂x∂t + C(x, t) ∂ 2 u ∂t 2 = Υ(u, ∂u ∂x , ∂u ∂t , x, t), (II.37)
where A, B and C are coefficients and Υ is an agglomeration of all operators of order less than 2. The PDEs are classified as

B 2 -4AC > 0, then equation (II.37) is hyperbolic (II.38) B 2 -4AC = 0, then equation (II.37) is parabolic (II.39) B 2 -4AC < 0, then equation (II.37) is elliptic. (II.40)
The behaviour of a PDE can be described using characteristics [START_REF] Zwillinger | Handbook of differential equations[END_REF]. These characteristics are ordinary differential equations that solve a given PDE along curves. These curves define the partial differential equation and demonstrate the general behaviour of the PDE during the solution phase.

II.2.1.1 Elliptic PDEs

Elliptic partial differential equations are steady state solutions of boundary value problems involving potentials [START_REF] Zwillinger | Handbook of differential equations[END_REF] [START_REF] Fletcher | Computational Techniques for Fluid Dynamics[END_REF], Fig. 2.9). The boundary values must be imposed on all four boundaries, so that the physical problem is well conditioned. The given potential problem will result in a smooth gradient field as depicted in figure II.1. A change in the boundary condition on any side is instantaneously felt by the exemplary point ①. This type of PDE is generally well-posed under some assumption on the regularity of the prescribed boundary values [START_REF] Fujita | Evolution problems[END_REF]. Its smoothing property disperses discontinuities, so that discontinuous solutions or numerical oscillation errors in general do not destabilise the solution [START_REF] Ascher | Numerical Methods for Evolutionary Differential Equations[END_REF]. The behaviour of parabolic partial differential equations is close to that of elliptic partial differential equations, in the sense that small perturbations immediately impact the solution everywhere in the domain. Parabolic PDEs describe transient diffusion problems. The heating equation (II.34) is of parabolic type, due to (II.39). A generic example for this class of partial differential equations is the one-dimensional homogeneous transient diffusion equation

∂u ∂t - ∂ 2 u ∂x 2 = 0. (II.42)
According to [START_REF] Hoffman | Numerical Methods for Engineers and Scientists[END_REF], this type of PDE has a primary temporal characteristic

t -t 0 = 0, (II.43)
as well as a secondary spatial characteristic

xx(t 0 ) = 0.

(II.44) A changing boundary condition for t 1 or t 2 is directly distributed to each point using the temporal characteristics ① and ②. Yet, its effect is inversely proportional to the distance of the boundary [START_REF] Hoffman | Numerical Methods for Engineers and Scientists[END_REF], such that the effect decreases with increasing distance. The characteristics ③ and ④ are spatial characteristics, which depend on the initial value u(t 0 ). E.g. the spatial point x 1 is influenced both at t 1 and t 2 by the boundary information, but also by its previous time increments. For the time increment t 2 it is, therefore, influenced by the temporal characteristic ① and by the spatial characteristic ④. This type of PDE is well-defined, if conforming boundary conditions are imposed and the complete boundary domain and initial values are prescribed, but numerical solution methods have to include the infinite information propagation speed [START_REF] Hoffman | Numerical Methods for Engineers and Scientists[END_REF].

II.2.1.3 Hyperbolic PDEs

According to [START_REF] Zwillinger | Handbook of differential equations[END_REF], a classical example for hyperbolic PDEs for which condition (II.38) holds, is the wave equation. The wave equation is defined as

∂ 2 u (∂t) 2 -c 2 ∇ • ∇u = 0, (II.45)
using a constant parameter c. This PDE resembles the parabolic heat equation, but features a second order derivative in time. In contrast to the above-mentioned elliptic and parabolic types of equations, it depends on the direction of the data flow, which is locally confined, so that local perturbations do not affect the complete computational domain. The general behaviour of this familiy of PDEs can be visualised using a special class of first order PDEs, which are the transient hyperbolic conservation laws [START_REF] Zwillinger | Handbook of differential equations[END_REF]. For a homogeneous one-dimensional setting, using the unknown quantity u, a general conservation equation can be defined as

∂u ∂t + v ∂u ∂x = 0, (II.46)
with v being a given velocity. The general form of this hyperbolic equation can be defined using the flux F (u) as

∂u ∂t + vF (u) = 0.
(II.47) [START_REF] Allen | Numerical analysis for applied science[END_REF] gives the characteristic curves for problem (II.46) using the initial time t 0 as

(x -x(t 0 )) -v(t -t 0 ) = 0. (II.48)
which are visualised in figure II.3.

Figure II.3: Characteristics for hyperbolic PDEs

Curve ① is the characteristic curve for the leftmost point of the domain defined at t 0 . The information of the inflow boundary condition flows along this line with velocity v. For this class of problems, the boundary condition must be consistent with the flow of information. It is visible that the information has only reached parts of the domain up to time t 2 . For this time step, the quantities of the redshaded domain only depend on the initial values x ≥ x m . It follows that a possible outflow boundary condition has no effect on the outcome of the computation. In effect, the simultaneous imposition of inflow and outflow boundary conditions can lead to inconsistencies, e.g. negative densities, due to numerical errors. A changing boundary condition will shift the characteristic in time, but it remains parallel to the initial curve, as can be seen for curve ② that starts from t 1 for the leftmost discretisation point. For the remaining spatial points the characteristic is shifted likewise for any given time step, as can be seen for curve ③ starting from x m . The characteristic curve ③ points in the direction of v, such that the left part of the domain is not influenced by the initial values presented in this part of the domain. In essence, for hyperbolic equations it is important that the boundary conditions are in accordance with the flow of information. The initial values do only have to be described, if the evolution for an initial quantity is needed. These initial quantities are not necessary to ensure the solvability of the partial differential equation. E.g. an open pipe of length L with a given inflow condition and a velocity v = L/t 1 will have filled in time t 1 irrespective of the initial quantities present in the pipe. The major drawback regarding hyperbolic PDEs is the dependency on information flow. A possible time discretisation for a numerical approximation is, therefore, not only restricted by the chosen discretisation method, but also by the physical problem itself. It is furthermore noticeable that a non-linear hyperbolic PDE can admit discontinuous solutions, even when smooth initial data is prescribed [START_REF] Ascher | Numerical Methods for Evolutionary Differential Equations[END_REF]. The weak formulation can admit non-physical solutions, which is why numerical solvers for hyperbolic conservation laws include filters that discard non-physical solutions. E.g. entropy filters can be used to discard quantities at jumps that would result in a loss of entropy. In essence, hyperbolic partial differential equations are difficult to handle, both based on a purely mathematical standpoint, regarding existence and uniqueness of solutions, as well as from a numerical point of view, due to stability and convergence issues. The general finite element method is, therefore, not commonly used in that context. For this class of PDEs additional work has to be done, in order to ensure solvability and stability. E.g. for advection type problems, the Streamline-Upwind-Petrov-Galerkin method [START_REF] Brooks | Streamline Upwind/Petrov-Galerkin Formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-Sokes Eequations[END_REF] is sometimes used, which creates an artificial diffusion term in the weak form to smooth discontinuities. Alternatively, the discontinuous Galerkin method is often used, since it allows naturally to introduce discontinuities into the discretisation as described in [START_REF] Cockburn | TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Conservation Laws II : General Framework[END_REF], [START_REF] Cockburn | The Local Discontinuous Galerkin Method for Time-Dependent Convection-Diffusion Systems[END_REF]] and [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF].

The E-H-J formulation, given in equations (II.22) and (II.23) is of hyperbolic type. This formulation can be rewritten in matrix formulation, such that the two curl-operators can be identified as flux (see equation ( 14) in [START_REF] Hesthaven | Nodal High-Order Methods on Unstructured Grids[END_REF]). The above-mentioned equation system has been extensively used in electromagnetic scattering applications. E.g. [START_REF] Hesthaven | Nodal High-Order Methods on Unstructured Grids[END_REF] solves the electromagnetic scattering problem using an explicit formulation using discontinuous Galerkin finite elements (cf. [START_REF] Hesthaven | Nodal discontinuous Galerkin methods: algorithms, analysis, and applications[END_REF]). For these problems, the objective is to model the travelling wave problem of electromagnetics. This is in contrast to induction heating, where it is assumed that the system of inductor and workpiece is in a quasi-steady state, so that the system itself behaves more like a diffusion problem, so that small perturbation immediately affect the solution values everywhere in the computational domain.

II.2.2 Solution methods and discretisation

The numerical solution of the eddy-current equations involves a discretisation step. The choice of discretisation depends on different factors. First, an assessment must be made regarding the computing resources that are available, since different methods possess different needs regarding computational time or memory capacity. Secondly, different methods behave differently with respect to numerical inconsistencies like oscillations, discretisation errors or problems with respect to floating point accuracy. In general, a discretisation method ought to be stable and consistent, so that results can accurately be produced with confidence. In addition, it is important to state, whether the underlying physical problem is fully defined using boundary values only, or whether volumetric source terms or non-homogeneous material behaviour necessitates a fully volumetric description.

In the context of electromagnetism two general approaches are widely used that are based on a domain based integration of the weighted residual of the PDE: the boundary element method (BEM) and the finite element method (FEM).

Let the general PDE be defined using an operator Υ with unknown solution vector u and right-hand side b, so that the residual r(u) can be given as

r(u) = Υ(u) -b = 0.
(II.49)

The common approaches are based on a weighted integration of the residual II.49 with the goal of minimising the error e h = uu h with respect to the discrete solution vector u h .

The BEM is based on the assumption that the underlying PDE can be transformed into a boundary integral equation using the fundamental form of the differential operator of the PDE. It is furthermore assumed that the enclosed computational domain consists of homogeneous materials, that it reacts linearly and that no inhomogeneous volumetric source excitations exist. An overview of this discretisation approach can be found in [START_REF] Sauter | Boundary Element Methods[END_REF]. The advantage of BEM is that only the boundary of the domain has to be approximated, which simplifies the application regarding complex domains. It also naturally permits relative movements of multiple domains. Applications are moving boundary problems in electromagnetics, e.g. [START_REF] Frangi | Magneto-mechanical simulations by a coupled fast multipole method-finite element method and multigrid solvers[END_REF] or [START_REF] Alotto | Three-dimensional eddy current analysis in unbounded domains by a DEM-BEM formulation[END_REF]. The disadvantage is that the surface integral approach necessitates a direct coupling of all elements, which results in a direct coupling of all degrees of freedom. For a practical application, this results in linear systems of equations that are difficult to treat numerically, due to their dense nature, even though their numbers of degrees of freedom are highly reduced.

The FEM is based on an elementwise approximation of the computational domain using discrete finite elements that, together, approximate the solution of the PDE. Applications for electromagnetics can be found in [START_REF] Kwon | Three-dimensional (3D) eddy current time-domain integral equation with Coulomb gauge condition[END_REF], [START_REF] Lu | Transient Eddy Current Calculation Using Magnetic Vector Potential Method Combined with Incomplete Tree Gauge[END_REF], [START_REF] Biro | On the Use of the Magnetic Vector Potential in the Finite Element Analysis of Three-Dimensional Eddy Currents[END_REF]] or [START_REF] Hömberg | A mathematical model for induction hardening including mechanical effects[END_REF]. Its discrete form leads to a global linear system of equations with a local coupling of the elements. Therefore, the discrete form is based on a sparse matrix system, which is simple to handle numerically. Furthermore, efficient solution methods can be employed, e.g. [START_REF] Hiptmair | Auxiliary Space Preconditioning for Edge Elements[END_REF] and [START_REF] Kolev | Parallel Auxiliary Sace AMG for H (curl) Problems. Technical report[END_REF].

A coupled approach, denoted as mixed FEM/BEM approach, can also be employed. In common electromagnetic applications, the embedding air-domain is discretised using a BEM approach, whereas the conductive domains are discretised using finite elements. The advantage of this coupling is that non-linear effects can be dealt with. Unfortunately, the resulting coupled linear system of equations is again difficult to treat numerically. Furthermore, the coupled linear system of equations is difficult to handle in parallel. Therefore, it can not efficiently be scaled, if the computational demand increases. Applications for this approach can be found in [START_REF] Hiptmair | Coupled boundary-element scheme for eddy-current computation[END_REF], [START_REF] Meddahi | An H-based FEM-BEM formulation for a time dependent eddy current problem[END_REF], [START_REF] Rodríguez | A FEM-BEM approach for electromagnetostatics and time-harmonic eddy-current problems[END_REF], [START_REF] Camaño | Analysis of a FEM -BEM model posed on the conducting domain for the time-dependent eddy current problem[END_REF] or [START_REF] Touzani | Mathematical Models for Eddy Currents and Magnetostatics[END_REF].

Figure II.4: Galerkin orthogonality principle

As stated above, both discretisation methods are applied in connection with an integration over the computational domain. Independently of the utilised method, the overall goal is to minimise the approximation error e h . In general, before integrating over the computational domain, the PDE is multiplied by a test function in vector form v, so that the result is in scalar form. For the BEM, a widely used method is the point collocation approach. Here, the test function v is chosen to be the Dirac delta function [START_REF] Gaul | Boundary Element Methods for Engineers and Scientists[END_REF], resulting in an evaluation of the underlying physical problem at these collocation points. Alternatively, the Galerkin approach, as described in [START_REF] Oden | An introduction to the mathematical theory of finite elements[END_REF] and [START_REF] Ciarlet | Handbook of numerical analysis[END_REF], can be employed, which chooses the test function space equal to the approximation space of the solution vector u. This approximation method is also the method most commonly associated with respect to finite element approaches.

As stated above, the Galerkin method is based on the assumption that both the approximation space U h of the solution vector u h , as well as the approximation space V h of the test function v h are equal. The Galerkin method possesses the advantage that the approximation error e h will be orthogonal in the approximation space, under the assumption that a solution u h exists, that it is uniquely defined and that the trial spaces are Hilbert spaces [START_REF] Girault | Finite Element Approximation of the Navier-Stokes Equations[END_REF]. This remarkable property is visualised in figure II.4. The orthogonality property indicates that the approximative solution u h ∈ U h must be the best possible approximation for the unknown vector u ∈ U in the full approximation space U h . It is, therefore, the optimal solution with the smallest error e h .

In addition, the weak form of the integral formulation can be used to reduce the order of the resulting formulation using integration by parts. In essence, Green's theorem [START_REF] Girault | Finite Element Approximation of the Navier-Stokes Equations[END_REF]] can be applied to the integrated form of the multiplied operator of the strong form of the PDE, in order to reformulate the problem into sums of volume and surface integrals. [START_REF] Oden | Finite elements: An introduction[END_REF] states that the solution is satisfied in a weak distributional sense. After certain assumptions on the test function space V, the surface integral formulations incorporating the higher order derivatives vanish, resulting in a reduction of the order of the PDE. In the following sections, this approach is applied to the strong form of the elliptic voltage potential problem (II.30), the parabolic heating equation (II.34) and the magnetic vector potential formulation (II.29).

II.2.2.1 Weak form of the voltage potential problem

The voltage potential problem, defined in equation (II.30) is a boundary value problem. Two different kinds of boundary conditions can be identified. A surface current J 0 • n can be prescribed on parts of the surface domain ∂Ω N , whereas a voltage potential can be prescribed on another part ∂Ω D .

∇ • σ∇Φ = 0,
in Ω (II.50)

σ∇Φ • n = J 0 • n, on ∂Ω N (II.51) Φ = Φ 0 , on ∂Ω D (II.52)
This boundary domain is based on the assumption that the complete surface domain is partitioned into these two subsets

∂Ω = ∂Ω N ∪ ∂Ω D , (II.53)
which do not overlap

∂Ω N ∩ ∂Ω D = ∅. (II.54)
The voltage potential problem is only uniquely defined, if boundary conditions are identified on the complete exterior domain and at least some part is identified as ∂Ω D . Equation (II.50) can be multiplied by a vector valued test function ṽ and integrated over the full domain Ω, with subsequent application of Green's theorem [START_REF] Girault | Finite Element Approximation of the Navier-Stokes Equations[END_REF], such that

 Ω ∇ • σ∇Φṽ dΩ = -  Ω σ∇Φ∇ṽ dΩ +  ∂Ω N (σ∇Φ • n)ṽ d∂Ω N +  ∂Ω D (σ∇Φ • n)ṽ d∂Ω D .
(II.55)

The normal derivatives of the voltage potential do not have to be explicitly defined on ∂Ω D , if the test function ṽ can be chosen, such that v = 0, on ∂Ω D .

(II.56)

In addition, by identifying the first surface integral of the right-hand side as the imposed current (II.51), equation (II.55) can be rewritten as

 Ω σ∇Φ∇ṽ dΩ    ã(Φ,ṽ) =  ∂Ω N (J 0 • n)ṽ d∂Ω N    b(ṽ)
.

(II.57)

The symmetric left-hand side of equation (II.57) is denoted by ã(Φ, ṽ). It is linear for argument ṽ for fixed Φ, as well as linear for argument Φ for fixed ṽ, such that it can be defined as a symmetric bilinear form. The right hand side b(v) can equally be identified as linear operator on argument ṽ.

It follows that the original PDE defined in equation (II.30) can be replaced by an equivalent boundary value problem using a weak formulation. It is only valid in a distributional sense [START_REF] Oden | Finite elements: An introduction[END_REF] . The advantage is that the PDE of second order with respect to the spatial coordinate, defined in equation (II.30) can be described by using the boundary value problem including only first order derivatives. The weak form of the voltage potential formulation is defined as: find Φ, such that ã(Φ, ṽ) = b(ṽ),

in

Ω σ∇Φ • n = J 0 • n, on ∂Ω N Φ = Φ 0 ,
-  Ω C ∇ • (k∇T )v dΩ C =  Ω C k∇T ∇v dΩ C -  ∂Ω CN (k∇T • n)v d∂Ω CN -  ∂Ω CD (k∇T • n)v d∂Ω CD .
(II.61) The convection flux (II.35) and the radiation loss (II.36) depend on the ambient temperature T A , as well as on the unknown temperature T , with the radiation loss being highly non-linear. According to [START_REF] Bay | A numerical model for induction heating processes coupling electromagnetism and thermomechanics[END_REF], it can be linearised using

q Fr = εσ s (T 4 -T 4 A ) ≈ εσ s (T t n-1 + T A )(T 2 t n-1 -T 2 A )    hr (T -T A ). (II.64)
The parameter h r is introduced explicitly, using the solution values of the previous time step T t n-1 . If the test function v is chosen, such that it vanishes on Ω CD and after inserting the heat flux (II.33) into equation (II.61), the bilinear form ā(T, v) of the weak form of equation (II.34) can be identified as

ā(T, v) =  Ω C ρCT v dΩ C +  Ω C k∇T ∇v dΩ C +  ∂Ω CN ((h + h r )T • n)v ∂Ω CN , (II.65)
whereas the linear form b(v) can be identified as

b(v) =  ∂Ω C Q em v dΩ C +  ∂Ω CN ((h + h r )T A • n)v d∂Ω CN .
(II.66)

The weak formulation of the initial value boundary problem with initial temperature T t 0 at time t 0 can, therefore, be stated as: find T , such that

ā(T, v) = b(v), in Ω C q F • n = (q Fr + q Fm ) • n, on ∂Ω CN T = T 0 , on ∂Ω CD T = T t 0 , for t = t 0 in Ω C .
(II.67) (II.68) (II.69) (II.70)

II.2.2.3 Weak form of the magnetic vector potential equation

The weak form of the magnetic vector potential equation can be written in a simplified manner, like the above-mentioned heat transfer equation, as well as the voltage potential equation. Integration by parts leads to Green's theorem including the curl operator (II.9) using vector valued u and v 

Ω (∇ × u)v dΩ =  Ω (∇ × u)(∇ × v) dΩ -  ∂Ω N (∇ × u × n)v dΩ N -  ∂Ω C (∇ × u × n)v dΩ C . (II.71)
The second term of equation (II.29) can, therefore, be rewritten in weak form after multiplication with a test function v and integration by parts over the domain Ω as

 Ω ∇ × 1 µ r ∇ × A dΩ =  Ω 1 µ r (∇ × A)(∇ × v) dΩ -  ∂Ω N ( 1 µ r ∇ × A × n)v dΩ N -  ∂Ω D ( 1 µ r ∇ × A × n)v dΩ D .
(II.72)

In the context of induction heating, the boundary surface ∂Ω N is mostly associated with symmetry conditions. In this work, it is assumed that the global domain Ω is chosen, such that it fully encloses the conducting surfaces with a size that is sufficient to assume that no magnetic flux exits or enters over this boundary surface, which can be identified as ∂Ω D . It is therefore assumed that the complete embedding boundary is equal to ∂Ω D , such that

∂Ω D = ∂Ω (II.73) ∂Ω N = ∅. (II.74)
It follows that, if the test function is chosen with ∇ × v = 0 on Ω D , the weak bilinear form of equation (II.29) can be written as

a(A, v) = ∂ ∂t  Ω σAv dΩ +  Ω 1 µ r (∇ × A)(∇ × v) dΩ, (II.75)
whereas the linear form can be written as

b(v) = -  Ω σ∇Φv dΩ (II.76)
under the assumption of a weak coupling between voltage potential Φ and magnetic vector potential A. The vanishing magnetic flux boundary condition can be prescribed on ∂Ω D as

B • n = 0 (II.77) (∇ × A) • n = 0, (II.78)
which can be fulfilled if

A × n = 0. (II.79)
The weak form of the initial value problem of the magnetic vector potential can thus be described using the initial value A t 0 at time t 0 as: find A, such that

a(A, v) = b(v), in Ω A × n = 0, on ∂Ω D A = A t 0 , for t = t 0 in Ω.
(II.80) (II.81) (II.82)

II.2.3 Finite element spaces

Finite element spaces are discrete representations of the vector spaces that are introduced in the weak Galerkin formulations (II.58), (II.67) and (II.80). These finite element vector functions are specified using the geometrical features of the chosen finite element discretisation. In practice, commonly used finite elements include triangles and quadrilateral elements in two dimensions, as well as tetrahedral, brick or wedge elements in three dimensions [START_REF] Wriggers | Nonlinear finite element methods[END_REF], [START_REF] Bathe | Finite Element Procedures[END_REF]. These discrete spaces are denoted with a subscript h, as

Φ ≈ Φ h (II.83) ṽ ≈ ṽh (II.84) T ≈ T h (II.85) v ≈ vh (II.86) A ≈ A h (II.87) v ≈ v h (II.88)
and replace the continuous representations in the respective weak formulations.

The representation choice depends on the modelling approach, as well as possible numerical and physical complications that might arise in a computational treatment. The interpolation function is based on a summation over the test functions with respect to the chosen geometrical features. In the following, an application will be given for triangular, as well as for tetrahedral elements, in order to detail the discretisation steps. A commonly used discretisation technique includes the Lagrangian finite elements [START_REF] Oden | An introduction to the mathematical theory of finite elements[END_REF], which are based on a nodal representation of the interpolating function. Let u be a generic representation for the different vector spaces (II.83)-(II.88). For a two-dimensional triangular element, featuring three independent nodes, a linear interpolation for the local elementwise form u e of u can be approximated as

u e (ξ, η) ≈ u eh (ξ, η) = 3  i=1 u i N i (ξ,
N i (ξ, η) = 0, if location (ξ, η) coincides with the local node index j ̸ = i. (II.91)
For the general linear triangle, defined in equation (II.89), the test functions can be chosen on a reference triangle, shown in figure II.5a, as such that conditions (II.90) and (II.91) are fulfilled. The advantage of the finite element method is that the geometry of each triangle element can be chosen in a flexible manner. Figure II.5b shows a real finite element, henceforth called the spatial configuration, which could be found in a finite element discretisation. Its orientation and form can be chosen arbitrarily. The reference finite element II.5a can be mapped onto the spatial configuration using a representation function T . This is the isoparametric concept [START_REF] Jin | The finite element method in electromagnetics[END_REF][START_REF] Wriggers | Nonlinear finite element methods[END_REF]. It is very advantageous for expressing the spatial derivatives of the local element function by applying the chain rule. Let ∂ ∂ξ and ∂ ∂η be the spatial derivatives with respect to the local coordinates of the reference finite element. A connection to the derivatives with respect to the global spatial coordinates x and y can be specified using the chain rule as

N 1 (ξ, η) = 1 -ξ -η (II.92) N 2 (ξ, η) = ξ (II.
 ∂ ∂ξ ∂ ∂η  =  ∂N 1 (ξ,η) ∂ξ ∂N 2 (ξ,η) ∂ξ ∂N 3 (ξ,η) ∂ξ ∂N 1 (ξ,η) ∂η ∂N 2 (ξ,η) ∂η ∂N 3 (ξ,η) ∂η  •    x 1 x 2 x 3       J •  ∂ ∂x ∂ ∂y  , (II.95)
by using the fact that the global coordinate x can be specified using the approach defined in equation (II.89) as

x(ξ, η) = 3  i=1 x i N i (ξ, η) =  N 1 (ξ, η) N 2 (ξ, η) N 3 (ξ, η)  •    x 1 x 2 x 3    .
(II.96)

The matrix J in equation (II.95) can be identified as the Jacobian transformation matrix. It can be derived for each distinct finite element to map between spatial and reference configuration.

The concept extends to finite elements of different shapes and dimensions. E.g. the linear tetrahedral finite element, which is based on four global coordinates

x =       x 1 x 2 x 3 x 4       ,
(II.97) can be specified using four local test functions, using

N 1 (ξ, η, ζ) = 1 -ξ -η -ζ (II.98) N 2 (ξ, η, ζ) = ξ (II.99) N 3 (ξ, η, ζ) = η (II.100) N 4 (ξ, η, ζ) = ζ. (II.101)
These linear nodal finite elements are discrete versions of the general Hilbert space H(Ω), involving square-integrable gradients, which can be defined as

H(Ω) = {u ∈ L 2 (Ω); ∇u ∈ L 2 (Ω)}.
(II.102)

Further test functions for higher polynomial orders and different element choices like quadrilateral elements or brick elements can be found in [START_REF] Wriggers | Nonlinear finite element methods[END_REF].

The choice of nodal values for describing the test functions is not unique. In fact, many more geometrical features can be used as reference for the description of the base functions. In electromagnetics it is found that the choice of Lagrangian nodal finite elements leads to spurious oscillations in the solution [START_REF] Jin | The finite element method in electromagnetics[END_REF], [START_REF] Bossavit | Computational electromagnetism: Variational formulations, complementarity, edge elements[END_REF]. This is due to the non-physical properties of the approximating Lagrangian finite element space with respect to electromagnetics. In fact, in electromagnetics only the tangential continuity must be enforced between interfaces, due to the use of the curl-operator (II.9) as stated in [START_REF] Jin | The finite element method in electromagnetics[END_REF], whereas in other physical problems like the heat distribution (II.67) or the voltage potential computation (II.58) the continuity condition is prescribed for the normal direction. A finite element space that introduces conforming basis function are the Nédélec/Whitney form finite elements [START_REF] Nedelec | Mixed Finite Elements in R3[END_REF], [START_REF] Bossavit | Computational electromagnetism: Variational formulations, complementarity, edge elements[END_REF]. For these finite elements the geometrical reference feature is chosen as the edge, i.e. the unknown, which is the magnetic vector potential A, is imposed on each edge degree of freedom. Since each edge k is uniquely defined by its border points i, j the edgewise test function Ψ k can be expressed using these nodes. For the simple linear tetrahedron, six local edge test functions can be defined as

Ψ k = Ψ ij = (N i ∇N j -N j ∇N i )L ij ,
(II.103) using the Lagrangian test functions of the border points i and j. The factor

L ij = ∥x i -x j ∥ (II.104)
is the edge length and has been included to normalise each test function, while making it dimensionless [START_REF] Jin | The finite element method in electromagnetics[END_REF]. The test functions are locally divergence free [START_REF] Nedelec | Mixed Finite Elements in R3[END_REF], [START_REF] Jin | The finite element method in electromagnetics[END_REF], so that

∇ • Ψ k = 0.
(II.105)

The resulting edge vector finite elements are a discrete version for a curl-conforming Hilbert space H(∇×, Ω) that can be defined as

H(∇×, Ω) = {u ∈ L 2 (Ω); ∇ × u ∈ L 2 (Ω)}. (II.106)
The weak formulations (II.58), (II.67) and (II.80) are integral formulas. Following [START_REF] Bathe | Finite Element Procedures[END_REF] and [START_REF] Wriggers | Nonlinear finite element methods[END_REF], these global integration formulas over the domain Ω are first replaced by a discrete version over T . Afterwards, the global integral is replaced by an assembly over the n e local finite element integrals, each having a unique volume T e , as (II.107)

It is advantageous to integrate each spatial element in the space of the respective reference element, such that the integral is transformed into the reference space.

For a three dimensional function f (x) the integral transformation is defined as

 Ωe f (x)dΩ e =  Ω ref f (ξ, η, ζ)det(J)dΩ ref , (II.108)
where Ω e is the volume of the spatial element, whereas Ω ref defines the volume of the reference element. The integration can be done numerically using Gauss quadrature rules [START_REF] Bathe | Finite Element Procedures[END_REF], such that equation (II.108) becomes a summation formula

 Ωe f (x)dΩ e ≈ np  p=1 f (ξ p , η p , ζ p )det(J(ξ p , η p , ζ p ))W p .
(II.109)

Here, W p is a weighting parameter of the quadrature rule for the reference domain, whereas ξ p ,η p and ζ p are the spatial coordinates in the reference space that are defined by the chosen quadrature rule, given the polynomial order p. [START_REF] Wriggers | Nonlinear finite element methods[END_REF] gives several parameters for commonly used finite elements. For a linear tetrahedral element it states that the integration is exact for a choice of p = 1, with ξ p = 1 4 , η = 1 4 and ζ = 1 4 and weighting factor W p = 1 6 .

II.2.4 Existence and uniqueness

It is mentioned in subsection II.2.2 that the Galerkin orthogonality results in the best possible approximations on the trial space U h for the real solution u ∈ U. An important concept in that regard is the existence and uniqueness of the solution vector. In fact, without uniqueness, the approximation is arbitrary and any statements regarding the residual error are essentially meaningless. For each of the above mentioned boundary value problems, namely the elliptic voltage potential equation (II.58) -(II.60), the heating equation (II.67)-(II.70) and the magnetic voltage potential formulation (II.80)-(II.82), existence and uniqueness can be defined under certain assumptions with respect to the bilinear operators and with regard to the regularity of the applied boundary and initial conditions.

The elliptic boundary value problem, defined in equations (II.58) -(II.60), is a classical example problem that has been studied in the literature. Many aspects with regard to existence, uniqueness and error control can be found in, amongst many others, [START_REF] Oden | An introduction to the mathematical theory of finite elements[END_REF] and [START_REF] Ciarlet | Handbook of numerical analysis[END_REF]. [START_REF] Oden | An introduction to the mathematical theory of finite elements[END_REF] states that the problem of existence and uniqueness, with respect to elliptic operators, can be reformulated using the Lax-Milgram-theorem [Lax and Milgram, 2005] as follows: If u ∈ U h and v ∈ V h are defined using inner product spaces, i.e. Hilbert spaces, then there exists a unique u for any arbitrary v, if the linear operator a(u, v) is bounded using some constant C 1 |a(u, v)| ≤ C 1 ∥u∥∥v∥, ∀u ∈ U and ∀v ∈ V (II.110)

and coercive using some constant C 2 a(u, u) ≥ C 2 ∥u∥∥u∥, ∀u ∈ U.

(II.111)

If the approximation space U h is then chosen to be equal to the trial space V h , the unique approximation u h will minimise the approximation error e h . For the context of electromagnetic induction [START_REF] Hömberg | A mathematical model for induction hardening including mechanical effects[END_REF] states that the above-mentioned Lax-Milgram theorem can be applied to the elliptic voltage potential problem. [START_REF] Hömberg | A mathematical model for induction hardening including mechanical effects[END_REF] states that the weak form (II.58) -(II.60) is uniquely solvable under the assumption that the inductor domain Ω I is closed, sufficiently regular and inside the computational domain Ω, any applied source currents are sufficiently smooth, square-integrable and regular and that the electrical conductivity is bounded inside

Ω I as 0 < σ < ∞. (II.112)
For the transient parabolic problem, existence and uniqueness are assured, if the bilinearform possesses the above-mentioned coercivity and boundedness for every fixed time step t in the bounded time interval t ∈ [0, t end ], with the additional assumption on sufficient regularity of the initial solution u 0 [START_REF] Lions | Remarques sur les équations différentielles opérationnelles[END_REF]. [START_REF] Hömberg | A mathematical model for induction hardening including mechanical effects[END_REF] provides conditions for the existence and uniqueness of the heating equation for a coupled thermomechanical phase transformation model that is weakly coupled to an electromagnetic model. [START_REF] Hömberg | A mathematical model for induction hardening including mechanical effects[END_REF] demonstrates existence and uniqueness of the temperature solution in weak form under the assumption that the thermal conductivity k, the density ρ and the specific heat capacity C are positive and bounded inside the workpiece domain

Ω C as 0 < k < ∞ (II.113) 0 < ρ < ∞ (II.114) 0 < C < ∞, (II.115) the heat transfer coefficient α is not negative α ≥ 0, (II.116)
the initial temperature T t 0 at time t = t 0 is at least Lebesgue integrable, i.e. T 0 ∈ L 1 (Ω C ) and the heating source contribution is square integrable

Q em ∈ L 2 (Ω C ).
The last regularity condition is an important aspect for the solution of the weak problem, defined in equation (II.67)-(II.70), since it provides for heating source contributions that are generally non-smooth due to a constant elementwise approximation.

The conditions for boundedness and coerciveness for the weak form of the parabolic magnetic vector potential problem (II.80)-(II.82), follows in the same way. Coercivity and boundedness is demonstrated by [START_REF] Hömberg | A mathematical model for induction hardening including mechanical effects[END_REF] under the assumption that in Ω the relative magnetic permeability is positive and greater or equal to one µ r ≥ 1 (II.117)

and the electrical conductivity complies with (II.112). It is furthermore assumed that the initial solution A t 0 is square-integrable, curl-conforming and divergencefree.

For the examples presented in this work, this is established by initialising the solution vector to zero for the first time step. The existence of a unique magnetic vector potential A for every subsequent time step t ∈ [0, t end ] is problematic, due to the fact that every solution can be perturbed by an additional term of a gradient space. Therefore, the magnetic vector potential admits infinitely many solutions. [START_REF] Hömberg | A mathematical model for induction hardening including mechanical effects[END_REF] remedies this problem by assuming that the source current is fully divergence free (cf. Homberg1999, Lemma 2.2). Under this assumption A is unique.

II.2.5 Stable time discretisation and weak coupling procedure

The transient magnetic vector potential problem, as well as the heat diffusion problem, defined in differential form in equations (II.29) and (II.34) have to be discretised in time. The choice of time discretisation is linked to the mathematical nature of the physical problem itself, e.g. the material behaviour or the shape of the imposed source, which can be either a source current for the magnetic vector potential or a heating source contribution for the heat diffusion equation. In addition, all time discretisation methods need to comply with the spatial discretisation of the continuous problems. Let h be a characteristic element size of an element of the suitable spatial discretisation that has been chosen for the weak forms (II.80), (II.67), let the parameter P represent the physical behaviour and let t be the time interval of the observation. The time discretisation methods can be broadly categorised based on the applicable time step size ∆t of the discretisation method and the computational cost of calculating every subsequent solution of the discretised linear system of equations. In general, the time step size can be defined as a function of the spatial discretisation, the mathematical nature of the partial differential equation and the time t ∆t = f (h, P, t).

(II.118)

Two general families of discretisation methods exist, which are single and multi-step methods. Multistep methods base the update of each time step on the solutions of previous computations, whereas one-step methods utilise only the solution at the previous time step to compute the updated solution. Each of these categories can be further divided based on the cost of the numerical application, e.g. storage cost and computational time or the complexity of the implementation itself.

The solution cost can in a first estimation be related to two different categories of time discretisation models, notably explicit and implicit methods. Explicit methods are characterised by the fact that a subsequent solution is a function of the previous time step only, whereas implicit methods relate the previous and subsequent time steps using a system of equations. In contrast to implicit methods, explicit methods are numerically cheap to apply, since they only necessitate an update of the solution variable of the previous time step. The disadvantage is in general that the size of the applicable time step ∆t is very small, when compared to implicit methods. Examples of multistep methods are the explicit Adams-Bashforth and the implicit Adams-Moulton scheme [START_REF] Allen | Numerical analysis for applied science[END_REF]. A widely used family of onestep methods are the Runge-Kutta (RK) methods that include the widely used implicit and explicit Euler methods and the trapezoidal rule. A general overview of these methods can be found in [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] or [START_REF] Ascher | Numerical Methods for Evolutionary Differential Equations[END_REF].

Both multistep and onestep methods have advantages and disadvantages. E.g. [START_REF] Allen | Numerical analysis for applied science[END_REF]] mentions that multistep methods are in general more economical for explicit schemes for comparable higher order onestep methods, since they are based on the solutions that are demanded by the user, whereas higher order onestep methods are based on intermediate time steps that are discarded after each global time step. A major advantage of onestep methods is given by [START_REF] Ascher | Numerical Methods for Evolutionary Differential Equations[END_REF] with respect to variable error control, which states that it is easier to implement variable step size methods using onestep RK schemes, due to the smaller overhead in the step size change with respect to the formulas of the multistep methods. In addition, higher order multistep methods necessitate an estimation of starting steps. These starting steps are often generated using RK methods or by using lower order multistep methods [START_REF] Hundsdorfer | Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations[END_REF].

Each of these methods creates approximations of the exact solution. Therefore, the approximation is by definition correct up to a certain approximation order. The application order relates the updated solution with the approximation error of the method. It is linked to the consistency of the numerical method. Let u t n+1 be the real solution and u t n+1 h the approximative solution. E.g. [Kraaijevanger, 1991] states that a solution method is consistent, if the error e(h) behaves like lim h→0 e(h) = 0, (II.119) with

∥u t n+1 -u t n+1 h
∥ ≤ e(h).

(II.120)

The order p ∈ N of the method can then be defined, if the approximation error of the numerical method behaves like e(h) = O(h p ).

(II.121)

An overview of the different approximation errors of common discretisation methods can be found in [Kraaijevanger, 1991], [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] or [START_REF] Ferracina | Strong stability of singly-diagonallyimplicit Runge-Kutta methods[END_REF]]. E.g. the application order of the standard explicit and implicit Euler schemes is p = 1, whereas the application order of the trapezoidal rule is p = 2.

Convergence is not a sufficient condition to ensure usability of a numerical approximation method. [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] states that stability analysis is vital to state the fitness of a numerical method. The behaviour of the numerical approximation with respect to numerical errors is important, to assess the applicability. [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] introduces the concept of stability using a model problem involving an unknown vector y and a linear operator M with the relation ∂y ∂t = My.

(II.122)

After introduction of the explicit Euler approximation scheme a recursive formula can be given, relating the n-th approximation vector y n to the initial value y 0 , using the identity matrix I as

y n = (I + ∆tM    R ) n y 0 , (II.123)
with R being identified as stability function. [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] Matrix M is chosen as the Jordan canonical form of M, i.e. M can be identified as a matrix in block form using the Jordan blocks J i , specified as

J i =        λ i 1 λ i 1 . . . . . . 1 λ i        .
(II.126) [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] notes that this transformation can be understood as a decoupling of the system of differential equations, since the other components of the Jordan block will depend on the solution, but will not influence its behaviour. It follows that (II.125) decouples for every general eigenvalue λ : λ i ∈ C, so that the stability function R can be written, using z = ∆tλ as

R(z) = 1 + z. (II.127)
For n → ∞ stability is assured, if the stability function (II.127) is bounded by 1. This is true for the condition

1 + z ≤ 1, (II.128)
which is visualised in figure II.6a. Analogously, [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] gives the stability function for the implicit Euler method as

R(z) = 1 1 -z . (II.129)
The resulting stability region is shown as shaded area visible in figure II.6b. [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] states that, when the stability function R(z) includes the full left half-plane of the space C of complex numbers, every approximation is bounded, if the exact function is bounded. This criterion is named A-stability and is according to [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] an essential property for numerical methods for approximating stiff PDEs. It can be seen that the explicit Euler scheme is not A-stable. This method thus possesses severe time restrictions due to the small stability region of admissible time step sizes. In fact, [START_REF] Ascher | Numerical Methods for Evolutionary Differential Equations[END_REF] and [START_REF] Iserles | A first course in the Numerical Analysis of Differential Equations[END_REF] state that no explicit RK method can be A-stable.

A further stability function can be stated for the trapezoidal rule, which is given by [START_REF] Leveque | Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems[END_REF] as

R(z) = 1 + 1 2 z 1 -1 2 z . (II.130)
It is thus defined as the full left half plane of the space of complex numbers C. The stability region is shown in figure II.6c. Additional stability regions for common RK methods and for the Adams-Bashforth and for the Adams-Moulton schemes are provided in [START_REF] Ascher | Numerical Methods for Evolutionary Differential Equations[END_REF] and [START_REF] Leveque | Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems[END_REF]. Following the definitions of [START_REF] Ehle | On Padé Approximations to the Exponential Function and A-Stable Methods for the Numerical Solution of Initial Value Problems[END_REF], oscillations of the solution are only damped, if the numerical method is A-stable and the stability function approaches 0, if the modulus z becomes large, i.e.

lim z→∞ ∥R(z)∥ = 0, (II.131)
which is then called L-stability. The trapezoidal rule is A-stable, but not L-stable, so that numerical oscillation errors are not damped out, which can lead to numerical instabilities.

It follows that for both the heat transfer equation, as well as for the magnetic vector potential problem, implicit L-stable methods, like the implicit Euler method should be employed. This ensures absolute stability, without restrictions on the applicable time step size. For higher approximation orders, [START_REF] Beck | Residual Based A Posteriori Error Estimators for Eddy Current Computation[END_REF] states that the semi-discrete magnetic vector potential problem can be discretised in time using any singly diagonally implicit Runge Kutta method (SDIRK) [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF], [START_REF] Ferracina | Strong stability of singly-diagonallyimplicit Runge-Kutta methods[END_REF].

Explicit schemes, should be avoided for these types of equations. In the context of electromagnetics, they are used in [START_REF] Hesthaven | Nodal High-Order Methods on Unstructured Grids[END_REF], which demonstrates the application of higher order explicit Runge Kutta schemes for Maxwell's equations in the hyperbolic E-H-J formulation, given in equations (II.22) and (II.23), for solving the scattering field problem of a travelling electromagnetic wave. The decoupling of the discontinuous elements enables a localised time stepping, so that the explicit nature of the time stepping scheme can be dealt with, without affecting the global time step size.

The underlying phsical problems of electromagnetism and heat transfer act on different time scales. A typical induction heat treatment procedure is accomplished in the range of 10 s-20 s, whereas the behaviour of the electromagnetic problem is linked to the frequency of the applied current load, which ranges roughly between 1 kHz and several MHz. A visualisation with respect to both time scales is provided in figure II.7. Generally, the discrete time steps for a general heat diffusion problem can be set to approximately 0.1 s to achieve a good balance between accuracy and computational time. Unfortunately, a direct numerical treatment of the coupled induction heating problem including electromagnetic effects is impossible. A remedy is proposed by [START_REF] Favennec | The ultraweak time coupling in nonlinear multiphysics modelling and related optimization problems[END_REF] which introduces an ultra-weak coupling procedure, replacing the strong interaction of electromagnetic and thermal effects using an averaging procedure on the heating power contribution. [START_REF] Favennec | The ultraweak time coupling in nonlinear multiphysics modelling and related optimization problems[END_REF] states that the relative changes of the heating power density of each electromagnetic period can be averaged for a characteristic time step of the heat diffusion computation. The general assumption is that the heat source contribution Q em can be replaced by a stable mean heat source contribution Qem that is averaged over a given electromagnetic period T * as

Qem = 1 T *  T * 0 σ∥ ∂A ∂t ∥ 2 . (II.132)
In order to account for non-linear effects during the electromagnetic computation, the heat source contribution is evaluated for several electromagnetic periods, until convergence to a user supplied convergence parameter ε c is reached. The electromagnetic computation is repeated for several full electromagnetic periods until the convergence criterion

Qem (T * n+1 ) -Qem (T * n ) Qem (T * n ) ≤ ε c (II.133)
is fulfilled. Here, T * n+1 denotes the current electromagnetic period, whereas T * n denotes the previous electromagnetic period. To account for non-linear effects during the electromagnetic computation, each period T * is dvided into n te time steps, depending on the level of accuracy, as

n te = 2 m , m = (4, 5, 6, . . . , n), for n ∈ N. (II.134)
After convergence has been established, it is then assumed that this stable mean heating power source Qem can be included into the linear form of the weak formulation (II.67) until the temperature induced phase changes, affecting the electromagnetic computation, i.e. conductivity σ and relative magnetic permeability µ r , necessitate a re-evaluation of the mean heating power density Qem .

II.3 Computation of conforming source currents

It is mentioned in section II.2.4 that the existence and uniqueness of the magnetic vector potential A strongly depends on the conformity of the imposed source current J. In fact, the existence and uniqueness theorem presented in [START_REF] Hömberg | Optimal shape design of inductor coils for surface hardening[END_REF] and [START_REF] Hömberg | A mathematical model for induction hardening including mechanical effects[END_REF] relies on the assumption that the source current forms a closed loop and is divergence free, so that

∇ • J = 0. (II.135)
It should be noted that the existence and uniqueness are based on the continuous representation of the weak form. For a discrete implementation, it can be assumed that discretisation errors and the general loss of accuracy due to the floating point representation lead to global inconsistencies with respect to the divergence of the magnetic vector potential, i.e. even a fully conforming system might possess small perturbations in the solution vector. With the use of Whitney form edge elements this error is minimised, since these elements are locally divergence free [START_REF] Nedelec | Mixed Finite Elements in R3[END_REF], [START_REF] Jin | The finite element method in electromagnetics[END_REF], even though the divergence free condition is not strongly enforced in a global manner (cf. [START_REF] Cockburn | Locally divergence-free discontinuous Galerkin methods for the Maxwell equations[END_REF]).

For the discrete problem, there exist two independently applicable methods to minimise the non-physical perturbation errors, due to gradients in the solution vector of the magnetic vector potential.

The first method is based on a graph theoretical approach for the discrete finite element mesh, using a tree-gauge. This tree-gauging method enforces the application on additional degrees of freedom in the interior of the domain, to ensure the divergence free condition of the magnetic field. It is based on the discrete representation of the inter-element facet flux using assumptions on the discrete representation of an element-wise conservation law for the magnetic vector potential A. In the following its influence on the solution and the convergence behaviour is presented for two test problems.

The second method enforces the correct application of the source current in a conforming manner, i.e. using discretisation and projection techniques, a possible divergence of the discrete source current is minimised. It is related to the initial computation of the source current on complex geometrical parts. Therefore, a possible solution technique for the introduction of physically conforming source currents is also described in this section. The conforming source current and its effect on the solution, as well as the convergence behaviour is demonstrated in detail for a three dimensional ring inductor test case. Finally, a combination of this technique, in connection with the tree-gauging method, is proposed.

The three dimensional test case consists of a cylindrical workpiece that is surrounded by a ring shaped inductor on which a voltage potential of Φ 0 = 10 V is prescribed. A geometrical representation is given in figure II.8a. The constant material data of the different domains, as well as the frequency of the applied current is given in table II

n periods f ρ in Ω I ρ in Ω c µ r in Ω I µ r in Ω C 16 1300 Hz 2 • 10 -8 Ωm 8 • 10 -7 Ωm 1 1
Table II.1: Material data and settings for the computation of the ring inductor test case

It is mentioned above that the source current has to be solenoidal, in order to conform to the mathematical model. In practice, this means that the source current has to form a closed loop. In the literature, this is mostly enforced by using simple circular shaped forms for the inductor geometry. E.g. [START_REF] Biro | On the Use of the Magnetic Vector Potential in the Finite Element Analysis of Three-Dimensional Eddy Currents[END_REF] proposes that the source current can be modelled as a source-contribution lying in the enclosing air-domain. [START_REF] Candeo | Multiphysics Modeling of Induction Hardening of Ring Gears for the Aerospace Industry[END_REF] demonstrates an induction hardening application for gearwheels. The source current is included using a singleturn induction coil, without taking into consideration the ring effect. [START_REF] Ren | Influence of the R.H.S. on the Convergence Behaviour of the Curl-Curl Equation[END_REF] shows that the form of the source current has a major impact on the solvability of the linear system of equations, demonstrating that small perturbations lead to non-physical source currents, which might lead to a loss of accuracy and loss of convergence properties with respect to the iterative solver. [START_REF] Ren | Influence of the R.H.S. on the Convergence Behaviour of the Curl-Curl Equation[END_REF] proposes to use the projection method presented in [START_REF] Albanese | Integral formulation for 3D eddy-current computation using edge elements[END_REF] to remedy this problem. In this work, these techniques are demonstrated for the example test case of the ring inductor.

II.3.1 Two cutting plane current computation technique

The computation of conforming source currents for inductors with complex geometrical features is not extensively dealt with in the literature. E.g. [START_REF] Ren | Influence of the R.H.S. on the Convergence Behaviour of the Curl-Curl Equation[END_REF] assumes that the source current follows a simple geometry that can be evaluated for every spatial coordinate, before the projection techniques are applied. Furthermore, it is mentioned in subsection II.2.2.1 that the weak form of the elliptic voltage potential problem is defined as boundary value problem. Therefore, boundary values must be defined on the free surface. Unfortunately, a closed geometrical shape can not support these surfaces directly, because it would invalidate the conformity of the solenoidality, i.e. a source current or a source voltage potential can not conformingly be prescribed inside of a finite element geometry. This is also noted by [START_REF] Touzani | Mathematical Models for Eddy Currents and Magnetostatics[END_REF], which states that "we are in presence of a paradox because a torus device does not allow for applying a current or voltage."

In the following, a simple algorithm is proposed to prescribe conforming and arbitrary source currents and voltage potentials that create conforming source currents on arbitrarily complex finite element discretisations. The key idea is the introduction of virtual cutting planes inside the discretised inductor geometry. A visualisation is provided in figures II.9a-II.9c. This approach is comparable to the approach of [START_REF] Touzani | Mathematical Models for Eddy Currents and Magnetostatics[END_REF], which also introduces a cutting plane to impose the data of the power source on a solenoidal inductor coil in a FEM/BEM context.

Let S 1 be a first virtual plane that is introduced, which cuts the inductor geometry, as shown in figure II.9a. It is bounded, so that it only cuts the geometry of the inductor one time. Based on this virtual plane, two element sets can be defined A surface source current or a surface voltage potential can now be readily applied on the closed inductor geometry using this discrete cutting surface ∂T -∩+ I by noting that a virtual copy of this set can be created with the exact same spatial coordinates and element connectivity, so that ∂T -∩+ Iv = ∂T -∩+ I . This virtual cutting plane creates additional degrees of freedom in the linear system of equations arising from a discretisation of the weak form, as defined in equations (II.58)-(II.60). For this first cutting plane, these additional degrees of freedom can be identified as reference Dirichlet voltage potential Φ 0v = 0 V. It is thus possible to prescribe an inner voltage potential Φ 0 ̸ = 0 V or a source current J • n on the cutting surface ∂T -∩+ I , while ensuring the solvability of the linear system of equations. The derived voltage potential can then be post-processed into a current that follows a closed loop.

An example of this proposed algorithm is given in figure II.10a for the Laplacian potential problem, using the continuous geometrical features

Ω I of T I , Ω -∩+ I of ∂T -∩+ I and Ω -∩+ Iv of ∂T -∩+ Iv , described as ∆U = 0, in Ω I (II.138) U = 1, on ∂Ω -∩+ I (II.139) U = 0, on ∂Ω -∩+ Iv .
(II.140)

Here, the example domain Ω I is defined as a disk of diameter 2, with a circular shaped hole with a diameter of 1. The vertical cutting plane S 1 is bounded to include only the top part of the domain, with a starting point at x = 0 and y = 0.5, defined so that its normal points in the x-direction. It serves as an example test case for a closed industrial inductor.

It can be seen in figure II.10a that the cutting plane accurately identified the set of virtual cutting nodes. The potential is resolved, leading to a smooth distribution of resulting quantities. The post-processed gradients are visible in figure II.10b. Initially, it seems that the resulting gradients follow a circular motion. Unfortunately, it can be seen in figure II.10c that the strong imposition of Dirichlet boundary conditions on the virtual cutting plane results in a discretisation error. The imposition of fixed values on the surface nodes of the virtual cutting plane results in gradients that do not follow the circular form of the example domain due to the linear approximation properties of the linear triangular finite elements, as is demonstrated in figures II.11a and II.11b. In fact, in the underlying finite element mesh, all Dirichlet boundary conditions are prescribed on the sides of the triangular finite elements, with most of the cutting surface defining a full triangle edge of the affected elements. It follows that for most of these surface elements, two of the three unknown degrees of freedom are prescribed. The value is equal on these nodes, such that a resulting gradient must point in the direction that is perpendicular to the imposed triangle edge, as shown in figure II.11b. The remaining degree of freedom affects the magnitude and sign of the gradient, but not its general direction.

A major advantage of elliptic PDEs is that they possess excellent smoothing properties, e.g. [START_REF] Fujita | Evolution problems[END_REF] mentions that even locally singular sources can be prescribed without destroying this smoothing property. The solution in figure II.10b indicates that this smoothing property results in correct gradients in most of the domain. A natural extension to the above-mentioned method consists therefore in the introduction of a second virtual cutting plane S 2 . For the presented example, the cutting plane has been chosen like S 1 but restricted to the bottom half of the geometry. The overall goal is to compute a potential problem on a closed domain, such that the initially prescribed potential is resolved. In assuming that the smoothing property of the elliptic PDE results in conforming values for the second cutting domain, it can be argued that it suffices to re-solve the original potential problem with these computed values. This is due to the fact that at this location, the gradient is following a correct circular path. Therefore, the nodal values of the existing second virtual cutting plane are correctly prescribed. It thus suffices to prescribe these values as Dirichlet conditions on the discrete nodes of the second cutting plane ∂T -∩+ I . The nodes on the resulting virtual cutting plane ∂T -∩+ Iv have to be chosen to conform to the initial potential difference, i.e. the nodal values of the virtual cutting plane have to be reduced by the value 1. A visualisation is provided in figure II.10d, where it can be seen that one cutting plane has prescribed values in the region of approximately 0.5, whereas the virtual plane possesses nodes with negative values of approximately -0.5, so that the initial potential difference remains. Figure II.10e visualises the post-processed gradient, which is essentially comparable to the initial gradient provided in figure II.10b. The key difference lies in the conformity of the circular gradients for the initial cutting surface. The resulting gradient is smooth and follows the circular direction. Figure II.10f provides the norm of the differences of the initial gradient and the recomputed gradient. the bottom half shows a good agreement in the gradients, whereas the top section, close to the first cutting plane, shows the strong non-conforming deviations of the first gradient.

It follows that this method can be used to compute conforming solenoidal source currents in arbitrarily complex closed inductor geometries. Figure II.12a presents the results of this two-cutting plane approach for an industrial inductor on which a voltage potential of Φ 0 = 25 V has been prescribed. The cutting planes used for prescribing the potentials are shown in figure II.12b. The current conforms to the inductor geometry, as can be seen in figure II.12c, which shows the streamlines of the post-processed electrical currents. The resulting currents can directly be used as source-contribution for the computation of the magnetic vector potential problem, described in equations (II.80)-(II.82). An example for a three-dimensional ring inductor test case is presented later in this section.

It must be noted that conforming internal surfaces could be generated by special meshing software, so that the initial computation directly generates conforming source currents. The drawback with this approach is that special meshing software is necessary to prescribe the initial conditions, whereas the described algorithm works for all types of finite element meshes. In addition, it must be noted that it can not always be known in advance, how the conforming internal surface should be placed, in order to derive solenoidal source currents, since the application depends on the complexity of the inductor geometry.

II.3.2 Tree gauging method

After the initial source current has been computed, the magnetic vector potential problem can be solved. As mentioned above, the first technique to improve the conformity of the resulting solution vector, while enforcing the divergence-free condition is the tree-gauging method.

The tree-gauging method is based on the discrete representation of the magnetic vector potential problem. The elemental description of the geometry, including inductor, workpiece and embedding air domain, is described using finite elements. These finite elements are described using vertices. Every element is defined by a local vertex adjacency matrix, which can be applied to the globally assembled linear system of equations. It follows that an undirected nodal adjacency graph can be found for the global finite element mesh that connects all nodal vertices. The connections can be identified as the edges of the finite elements. [START_REF] Albanese | Integral formulation for 3D eddy-current computation using edge elements[END_REF] states that the finite element description can not support divergences in the solution, if some nodal spanning tree is added to the imposed Dirichlet conditions. Early applications of this spanning tree technique were used to reduce the non-definiteness of the linear system of equations. In fact, [START_REF] Manges | A generalized tree-cotree gauge for magnetic field computation[END_REF] states that the tree-gauge ensures that the linear system of equations has a null space that is restricted to zero.

Figure II.13a shows a spanning tree for the ring inductor test case that has been created with the algorithm provided in [START_REF] Prim | Shortest connection networks and some generalizations[END_REF]. It is visible that the spanning tree connects all vertices of the finite element mesh, while avoiding cycles. The edges that are marked by the spanning tree algorithm are prescribed to a fixed value during the computation of the magnetic vector potential A, in order to avoid the creation of divergences. The remaining edges make up the co-tree, which is visible in figure II.13b. For an edge finite element description, the co-tree is thus the set of remaining unknown degrees of freedom. Figure II.13c visualises the cylindrical part of the workpiece domain that can be found inside of the embedding air domain. It is again visible that cycles are avoided, while every vertex node is connected at least once. The tree-gauging method has been used by many authors, e.g. [START_REF] Lu | Transient Eddy Current Calculation Using Magnetic Vector Potential Method Combined with Incomplete Tree Gauge[END_REF], [START_REF] Biro | On the Use of the Magnetic Vector Potential in the Nodal and Edge Finite Element Analysis of 3D Magnetostatic Problems[END_REF] and [START_REF] Ménach | Determination and Utilization of the Source Field in 3D Magnetostatic Problems[END_REF]]. [START_REF] Lu | Transient Eddy Current Calculation Using Magnetic Vector Potential Method Combined with Incomplete Tree Gauge[END_REF] solves the eddy current problem using an edge finite element approach in connection with an incomplete tree gauging method. The incomplete gauging is introduced, to ameliorate the conditioning of the linear system of equations, in order to increase the convergence rate of the iterative solver. [START_REF] Lu | Transient Eddy Current Calculation Using Magnetic Vector Potential Method Combined with Incomplete Tree Gauge[END_REF] reports that the incomplete tree-gauging method converges better than the full tree-gauging approach, but that the convergence fails at some point. It is further mentioned that the results are non-unique. [START_REF] Biro | On the Use of the Magnetic Vector Potential in the Nodal and Edge Finite Element Analysis of 3D Magnetostatic Problems[END_REF] states that the numerical stability depends on the choice of tree that is used to gauge the linear system of equations. Since the choice of tree is quite arbitrary, it is mentioned that there seems to be no ideal choice. [START_REF] Biro | On the Use of the Magnetic Vector Potential in the Nodal and Edge Finite Element Analysis of 3D Magnetostatic Problems[END_REF] concludes that a simple edge finite element approach without gauging condition should be chosen over a tree-gauged version, since the accuracy is comparable, but convergence is faster. [Le [START_REF] Ménach | Determination and Utilization of the Source Field in 3D Magnetostatic Problems[END_REF]] imposes a conforming source current by utilising a continuity approach based on the discrete flux over the element facets, so that the current is divergence free. A tree gauging condition is included in this pre-processing step. It is mentioned that the convergence is faster than for a gauged linear system of equations. Yet, the results indicate that the convergence can only be achieved up to a certain residual error, like for the test cases presented in [START_REF] Lu | Transient Eddy Current Calculation Using Magnetic Vector Potential Method Combined with Incomplete Tree Gauge[END_REF]. The same approach is used by [START_REF] Dular | A Generalized Source Magnetic Field Calculation Methof for Inductors of any Shape[END_REF], but it it is not mentioned, what kind of effects it has on the convergence rate.

The impact on the convergence behaviour for the ring inductor test case is demonstrated in figures II.14a and II.14b. Both figures show the convergence behaviour of a preconditioned GMRES solver for the first time step of the electromagnetic computation with respect to the reduction of the relative residual errors, which are scaled to enable a direct comparison. In each case, the source current is prescribed using the above-mentioned two-cutting plane technique. The blue lines show the convergence of the ungauged linear system of equation, whereas the red lines indicate the behaviour of the gauged linear system of equations.

The preconditioner used for the test case presented in figure II.14a is an algebraic multigrid preconditioner using a Cleary-Luby-Jones-Plassmann (CLJP) coarsening technique [START_REF] Henson | BoomerAMG : A parallel algebraic multigrid solver and preconditioner[END_REF]. It is visible that the initial ungauged version of the solver converges rapidly, with a huge reduction of the relative residual error to 10 -7 in less than 20 GMRES iterations. The introduction of the treegauging results in a large decrease of the convergence rate. The convergence curve is considerably flatter, leading to a comparable reduction of the relative residual error for the ungauged version in approximately 200 iterations. The impact of the tree-gauging technique is, therefore, a reduction for the convergence rate by a factor of 10.

Figure II.14b visualises the convergence behaviour using a Jacobian preconditioner. The convergence rate is equally slow for both the gauged, as well as for the ungauged linear system of equations, as compared to the previous algebraic multigrid preconditioner. Even though the tree-gauged version initially converges better, it can be seen that the reduction of the residual error to 10 -4 takes considerably longer, taking into account that the plotted graph is visualised using a double-logarithmic scale. The final residual error that has been computed after 2000 GMRES iterations is slightly lower for the tree-gauged version, even though it remains still higher than the residual error of the algebraic multigrid version. The results indicate that the convergence behaviour can not be determined in advance, since the same tree-gauge leads to slightly different convergence behaviour for two different preconditioners. This is in accordance with [START_REF] Biro | On the Use of the Magnetic Vector Potential in the Nodal and Edge Finite Element Analysis of 3D Magnetostatic Problems[END_REF], which mentions that the arbitrary choice of the tree leads to inconsistent applications. A remedy is proposed in [START_REF] Hiptmair | Multilevel Gauging for Edge Elements[END_REF] that proposes an algorithm for computing almost optimal tree-gauges. Unfortunately, it is based on a hierarchical set of nested finite elements and is, therefore, difficult to apply in practice. Therefore, the results indicate that the global tree-gauging method should be avoided, due to the problematic impact on the convergence rate.

II.3.3 Current potential formulation

The second technique for reducing non-conforming divergences in the solution vector for the magnetic vector potential is the source current projection method, presented in [START_REF] Albanese | Integral formulation for 3D eddy-current computation using edge elements[END_REF]. [START_REF] Albanese | Integral formulation for 3D eddy-current computation using edge elements[END_REF] shows that the source current J can be represented by a vector potential formulation, using the source current vector potential H, so that ∇ × H = J.

(II.141)

The advantage of this formulation is that an application of the curl operator on the vector potential H automatically filters the solution, like the post-processing of the magnetic flux B from the magnetic vector potential A. [START_REF] Albanese | Integral formulation for 3D eddy-current computation using edge elements[END_REF] states that Nédélec vector type elements can be used to describe this vector potential, since the tangential inter-element continuity ensures the conformity of the normal components of the post-processed source current J. The addition of the large null space in the linear system of equations leads to non-uniqueness of this vector potential. Therefore, a tree-gauging can be employed to ensure the uniqueness of the vector potential. [START_REF] Ren | Influence of the R.H.S. on the Convergence Behaviour of the Curl-Curl Equation[END_REF] demonstrates the effects of this method on the convergence behaviour of an iterative solver for a symmetrical coil on which a solenoidal current is applied. In the beginning, [START_REF] Ren | Influence of the R.H.S. on the Convergence Behaviour of the Curl-Curl Equation[END_REF] shows that the form of the applied source current has a strong impact on the convergence behaviour. The initial test includes a source current that is only approximatively sinusoidal, with an application of the current in straight segments. The resulting non-physical source current distribution leads to convergence problems, so that the residual error can not be reduced further after a certain number of iterations. The second test includes a fully circular source current, applied to the coil geometry. The convergence is better, but the residual can again only be reduced up to a certain constant. Afterwards, it is shown that an application of the equation (II.141) leads to a strong improvement in the convergence rate. The residual error can be reduced, with increasing number of solver iterations. [START_REF] Ren | Influence of the R.H.S. on the Convergence Behaviour of the Curl-Curl Equation[END_REF] demonstrates that the vector potential form H of the source current can be computed using a weak formulation, by projecting the source current J on the curl of the test function space v ∈ V ⊂ H(curl, Ω), formulating the problem as: there exists

H ∈ U ⊂ H(curl, Ω), such that  Ω (∇ × v)(∇ × H) dΩ =  Ω (∇ × v)J dΩ, ∀v ∈ V. (II.
(v) =  Ω v(∇ × H) dΩ.
(II.143) [START_REF] Ren | Influence of the R.H.S. on the Convergence Behaviour of the Curl-Curl Equation[END_REF] reports an improvement of the convergence rate if the projection of equation (II.143) is reformulated using integration by parts,

b(v) =  Ω (∇ × v)H dΩ.
(II.144)

II.3.4 Benchmark application for a ring inductor test case

In the following the impact on the convergence rate is demonstrated for different combinations of the above-mentioned techniques for the ring inductor test case that is presented in figure II.8a. The convergence history is visualised for the computation of the magnetic vector potential using different forms of the right-hand source current. It is mentioned above that the one-cutting plane technique leads to non-conforming source currents. Therefore, in the following, the impact on the convergence, as well as the solution of the magnetic vector potential with respect to conforming and nonconforming source currents is visualised.

First, the source current is computed using the elliptic voltage potential formulation (II.58)-(II.60) using the virtual cutting plane technique. Figures II.15a and II.15c show the potential computation for one and two cutting planes. Figure II.15a shows the voltage potential with a voltage potential difference of 10 V,which is applied using a cutting plane in the top region of the ring inductor using a normal vector pointing in the x-direction. The two cutting plane test case is visible in figure II.15c with the second cutting plane being applied to the left region of the inductor, with a normal pointing in the y-direction. The voltage difference on the cutting nodes and the virtual cutting nodes equals the initial voltage potential difference of 10 V.

The post-processed source currents are shown in figure II.15b for the one cutting plane test case and in figure II.15d for the two-cutting plane test case. The current distribution and direction is visualised using vectors. It can be seen that the single cutting plane technique leads to inconsistencies with respect to the source current near the first cutting plane. The disturbances are non-physical and due to the above-mentioned discretisation error, resulting from the strong enforcement of the Dirichlet voltage conditions on the finite element boundaries. In contrast, the source current for the two cutting surface technique, shown in figure II.15d, follows the geometry of the ring inductor. There are no inaccurate approximation errors near the first nor the second cutting plane. Therefore, it can be argued that the source current, visible in figure II.15d, ought to be used in the computation of the magnetic vector potential A, whereas the current in figure II.15c should be avoided. In the following, the effects of the use of both currents are demonstrated with regard to the solution vector A, as well as the general convergence behaviour. In addition, the effects of different right hand side configurations on the resulting post-processed magnetic flux B are shown. The primary goal is to visualise the effects of non-physical source currents both on the vector potential A, as well as on the magnetic flux B. The use of the above-mentioned source current projection technique and subsequent use of the projected vector potential for the right-hand side of the magnetic vector potential equation is highlighted. Figures II.17 the results for a projected vector potential current H that has been computed after the global linear system of equations, resulting from a discretisation of equation (II.142), is gauged using a nodal spanning tree in order to ensure uniqueness of the potential with subsequent projection using equation (II.144). The figures II.17 and II.18 show the magnetic vector potential for the first time step of the magnetic vector potential computation for two different configurations of linear solvers. Figure II.17 shows the results that are obtained after 2000 Jacobian preconditioned GMRES solver iterations, whereas figure II.18 gives the results for an algebraic multigrid preconditioned GMRES solver using a CLJP coarsening technique. The respective convergence behaviour for each test case is given in figures II.16a-II.16d. Figure II.16a shows the convergence behaviour for the algebraic multigrid preconditioned solver, with a close up given in figure II.16b. The convergence behaviour of the Jacobian preconditioned solver is given in figure II.16c with an overview of the last iterations given in figure II.16d. It is assumed that the tree-gauged, two cutting surface technique leads to the most accurate results. It can furthermore be seen that the algebraic multigrid preconditioned solver converges rapidly, with a large reduction of the residual error, whereas the Jacobian fails to converge to the same residual error in the maximum number of iteration. Therefore, all results are scaled by the maximum value that is computed for the configuration using an algebraic multigrid preconditioner together with a source current contribution computed by a tree gauged vector potential formulation. This enables a direct visual comparison of the magnetic vector potential A, as well as of the post-processed magnetic flux B.

Both figures II.17 and II.18 give qualitatively comparable results. For both figures the color scale is chosen to reflect outliers in the computation. The scale is chosen, so that the magnetic vector potential, which is scaled by the maximum value of the tree-gauged two cutting plane projection using H is directly comparable. The scale is slightly enlarged to include values that are up to 10% larger than the reference value. Every quantity exceeding this limit is visualised in black. For both figures it is visible that the one cutting plane technique with direct application of the computed source current, shown in subfigure (a), leads to strong inconsistencies. The values close to the cutting plane strongly exceed the values to be found for the two cutting plane technique in subfigures (d)-(f). A small location near the cutting plane exceeds the comparison value by less than 10%, whereas a large part exceeds this value. The introduction of the projection method, using equation (II.144) effectively filters the most extrem values, as can be seen in subfigure (b). Yet, the computed value is lower than the computed reference value in subfigure (f). An additional tree-gauging, in order to ensure uniqueness of the computed vector potential H for the source current does not qualitatively alter the results. For both figures II.17 and II.18 it can be seen that the non-physical source current leads to deviations in the magnetic vector potential, even in the workpiece. The computed field is slightly excentric, deviating in the positive y-and negative x-direction.

Following from the magnetic vector potential A, the magnetic flux B can be post-processed. From the observation that the magnetic vector potential in figures II.17 and II.18 are comparable, it follows that it suffices to post-process the results of the algebraic multigrid preconditioned solution technique, in order to visualise the magnetic flux density B. Table II.2 gives an overview of the maximal differences for each computation. The reference value is again the two cutting surface test case using a tree-gauged H projection technique in connection with an algebraic multigrid preconditioned GMRES solver. For both the magnetic vector potential A, as well as the magnetic flux density B, the maximum relative difference of the maximal values is computed. For the one cutting plane technique it is visible that the maximal value exceeds the comparison by more than 400%, whereas the H projection technique effectively filters these maximum value, so that the maximum value is very close to the comparison value. The results for the magnetic flux B demonstrate again the influence of non-physical source currents. The initial result for the one cutting plane technique leads to a maximum value that is almost 700% as large as the maximal value of the comparison test case. Interestingly, the projected current, which leads to comparable magnetic vector source fields, leads to strong inconsistencies with respect to the magnetic flux. The absolute difference for both projected one cutting plane test cases exceeds the maximum value by more than 300%. In comparison, all two cutting plane test cases lead to comparable computational results, both for the magnetic vector potential, as well as for the post-processed magnetic field. for the magnetic vector potential and the magnetic field for the first time step of the electromagnetic computation with respect to a preconditioned GMRES solver using a Jacobian and an algebraic multigrid preconditioner.

∥ max( Ã)-max(A) max( Ã) ∥ • 100% ∥ max( B)-max(B) max( B)
The impact of the above-mentioned current contribution techniques on the convergence behaviour is demonstrated in figure II.16. It can be noted that the algebraic multigrid preconditioned solver, shown in figures II.16a and II.16b, converges rapidly, in much less iterations than the Jacobian preconditioned solver. In fact, the Jacobian preconditioned solver failed to reach a comparable relative residual error in the maximum number of iterations. It can further be noted that even though the non-physical test cases with non-conforming source currents converged rapidly, the resulting magnetic vector potential is essentially meaningless and not useful for a further treatment. The Jacobian preconditioned solver, shown in figures II.16c and II.16d, initially converges very slowly. The convergence rate increases after approximately 400 iterations and flattens considerably after 1000 iterations.

The algebraic multigrid method converges essentially in between 13 and 15 iterations. It can be argued that, therefore, the order of convergence is not distinguishable. In contrast, the double logarithmic scale of the Jacobian preconditioned test case enables a direct comparison of each source contribution technique. The shows that the remaining residual error is greatest for the non-physical contribution of the one cutting plane technique. The relative residual error is smaller for all the currents contributed using the two cutting plane technique. Naturally, the projected vector potential approach using the vector potential H leads to the greatest reduction of the residual error, which can be reduced further by introducing a gauging technique.

The results indicate the validity of the two cutting plane approach for the source current calculation of complex geometrical parts. Furthermore, it shows that the projection technique of [START_REF] Albanese | Integral formulation for 3D eddy-current computation using edge elements[END_REF] together with the weak formulation of [START_REF] Ren | Influence of the R.H.S. on the Convergence Behaviour of the Curl-Curl Equation[END_REF] leads to conforming source contributions. In contrast to the tree-gauging restriction of the global linear system for computing the magnetic vector potential A the tree-gauging is limited to the right-hand side contribution. Therefore, there is no adverse affect on the convergence of the linear system of equations nor with respect to the quality of the result. A reduction of the degrees of freedom is, therefore, avoided. The application can be rapidly done, since the projection technique can be based on an existing edge finite element discretisation. The resulting weak form, discretised using linear Nédélec vector finite elements, can be resolved using standard solution techniques. The examples indicate that a treegauging technique of the source current reduces the residual error. Yet, the results also indicate that for a practical application it suffices to project a conforming source current on a vector potential H, in order to ensure the conformity of the resolved magnetic vector potential A.

III Efficient linear solvers for the associated electromagnetic problem

Français:

La solution répétée de la formulation en potentiel vecteur magnétique nécessite l'application des méthodes de solution efficaces.

English:

The repeated solution of the electromagnetic magnetic vector potential formulation necessitates efficient solution methods. The following chapter gives an overview of some choices regarding classical Krylov subspace solution methods. The preconditioning of the linear system of equations is an important aspect of the solution phase. The preconditioning operator must be efficient, scalable and easy to setup. An overview of some classical preconditioning techniques is given. It is followed by the introduction of the auxiliary space preconditioning technique. This method is based on an L2-orthogonal Helmholtz-type splitting of the discrete curl-conforming finite element space. The transfer operators are easy to setup and apply. The resulting preconditioners give rise to elliptic bilinear forms and can be used to create efficient solution methods based on the multigrid method. Hence, the algebraic multigrid method can be used, so that the application is not dependent on special meshing software. The operators are based on the nodal description, so that standard algebraic multigrid solution techniques can be employed. An important aspect in that regard is the coarsening method of the multigrid operator. An overview of common methods is provided. The chapter demonstrates the application of efficient solution techniques for a problem involving a crankshaft in a global edge finite element configuration. Several benchmark results for different solver configurations are provided. They demonstrate the efficiency of the auxiliary subspace preconditioning technique in connection with a monotonically converging GMRES method.

III.1 Requirements for a linear solver for induction heat treatment

In this work the numerical model of induction heat treatment consists of two weakly coupled parabolic equations that are discretised in space using finite elements in connection with the use of implicit methods in time. Therefore, there exists a need for the repeated solution of the assembled linear systems of equations. The heat diffusion equation is defined for the conducting workpiece domain Ω C . It is discretised using nodal Lagrangian finite elements. Classical solution techniques, like the preconditioned conjugated gradient method (CG) in connection with a simple diagonal scaling Jacobian preconditioner suffice to achieve rapid convergence. The discretised linear system of equations is positive definite and well conditioned. Therefore, it is easily invertible, so that a classical LU-decomposition technique (see, e.g. [START_REF] Kelley | Iterative Methods for Linear and Nonlinear Equations[END_REF] or [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]) can be employed, if the memory space is sufficiently large. Alternatively, a Cholesky-decomposition technique (see, e.g. [START_REF] Kelley | Iterative Methods for Linear and Nonlinear Equations[END_REF] or [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]) can be employed, since the linear system of equations is symmetric. An incomplete LU-decomposition technique can naturally be employed, if the application of the full LU-decomposition technique would exceed the computational resources. Furthermore, due to the strong ellipticity of the heat diffusion problem, it is very advantageous to implement multigrid solvers (see, e.g. [START_REF] Brandt | Guide to multigrid development[END_REF]) or algebraic multigrid solvers that are based on the global adjacency information of the nodal element connectivity (see, e.g. [START_REF] Henson | BoomerAMG : A parallel algebraic multigrid solver and preconditioner[END_REF]). Implementations are often given for the standard Poisson benchmark model, but results are comparable to the transient parabolic heating equation. Some benchmarks problems are presented in [START_REF] Henson | BoomerAMG : A parallel algebraic multigrid solver and preconditioner[END_REF] and [START_REF] Stüben | A review of algebraic multigrid[END_REF].

In contrast, the global linear system of equations of the magnetic vector potential problem (II.80) -(II.82) is generally more difficult to treat. The static form of the magnetic vector potential problem is not as easily solvable as the static form of the heat diffusion equation, i.e. the Poisson equation. The reason for this is that the discrete version of the curl-curl operator, found in the bilinear form of equation (II.80), possesses a large kernel [START_REF] Manges | A generalized tree-cotree gauge for magnetic field computation[END_REF] that has to be addressed during the solution phase. This results from the fact that any gradient potential can be added to the resulting magnetic vector potential, without invalidating equation (II.25). In fact, [START_REF] Hiptmair | Nodal Auxiliary Space Preconiditioning in H(curl) and H(div) Spaces[END_REF] states that the curl-operator results in a linear system of equations with an infinite null space, since the full space of gradient functions has to be included. A classical remedy to this problem is the introduction of a gauging condition

∇ • A = 0.
(III.1) [START_REF] Albanese | Integral formulation for 3D eddy-current computation using edge elements[END_REF] introduces this gauging condition directly, based on a graph theoretical approach applied in connection with the shape functions of the finite elements. [START_REF] Albanese | Integral formulation for 3D eddy-current computation using edge elements[END_REF] proposes that the degrees of freedom should only be described on the edges that form a co-tree of the discrete computational domain, since every co-tree edge closes a unique loop, so that Stoke's theorem results in a uniquely defined potential field. The resulting gauged linear system of equations is invertible. [START_REF] Manges | A generalized tree-cotree gauge for magnetic field computation[END_REF] shows that this partitioning of the finite element space into tree and co-tree edges results in two orthogonal subspaces. Unfortunately, the tree-gauging is not uniquely defined, since the defining tree graph can be almost arbitrarily defined [START_REF] Biro | On the Use of the Magnetic Vector Potential in the Nodal and Edge Finite Element Analysis of 3D Magnetostatic Problems[END_REF]. Therefore, its application is not consistent, which can lead to severe convergence problems in connection with iterative solvers, as demonstrated in [START_REF] Ren | Influence of the R.H.S. on the Convergence Behaviour of the Curl-Curl Equation[END_REF] and [START_REF] Golias | Magnetostatics with Edge Elements: A Numerical Investigation in the Choice of the Tree[END_REF]. A method for the optimal selection of the tree is proposed in [START_REF] Golias | Magnetostatics with Edge Elements: A Numerical Investigation in the Choice of the Tree[END_REF]. It shows that an optimised spanning tree can be created, based on an elimination of the edges that lie in the principal directions of the magnetic vector potential A. Initially, a first approximative source potential distribution is generated, based on the source current J, which is projected onto the edge degrees of freedom. Afterwards, the spanning tree is generated, neglecting the trees with large imposed values. The resulting co-tree is optimised, so that the principal direction of the magnetic vector potential is covered by the unknown degrees of freedom. [START_REF] Golias | Magnetostatics with Edge Elements: A Numerical Investigation in the Choice of the Tree[END_REF] demonstrates the effects on sample applications, noting that the use of an arbitrary spanning tree reduces the numerical accuracy and leads to convergence issues with respect to an incomplete Cholesky type preconditioned CG method. The optimal source tree is found to increase the convergence rate and numerical accuracy, even though the ungauged version is still slightly more accurate and converges faster than the aforementioned techniques. [START_REF] Dlotko | Efficient generalized source field computation for h-oriented magnetostatic formulations[END_REF] and [START_REF] Biro | On the Use of the Magnetic Vector Potential in the Nodal and Edge Finite Element Analysis of 3D Magnetostatic Problems[END_REF] state this efficiency problem as the main reason for using an ungauged formulation. [START_REF] Lu | Transient Eddy Current Calculation Using Magnetic Vector Potential Method Combined with Incomplete Tree Gauge[END_REF] states that in practice it suffices to establish convergence up to a certain numerical error, supposing that the resulting magnetic vector potential is resolved accurately enough to be usable in further applications. The argumentation of [START_REF] Beck | Residual Based A Posteriori Error Estimators for Eddy Current Computation[END_REF] follows in the same way stating that the quantity of interest is the post-processed curl of the computed vector potential, which in this work is the magnetic field B, i.e.

B = ∇ × A.
It is stated that the magnetic flux density will remain unique, even if the computed magnetic vector potential is non-unique. Yet, the induction heating formulation in this work is identified based on the direct application of the eddy current equation in vector potential form, using the transient equation (II.26), so that the induced eddy currents inside Ω C are given by

J = -σ ∂A ∂t .
Therefore, it is important to compute a correct magnetic vector potential. In that regard, it is established in subsection II.2.4 and II.3 that the conformity of the imposed right-hand side is an important step to ensure the uniqueness of the vector potential. It is furthermore shown that the direct tree-gauging, as defined in [START_REF] Albanese | Integral formulation for 3D eddy-current computation using edge elements[END_REF] or [START_REF] Manges | A generalized tree-cotree gauge for magnetic field computation[END_REF] leads to strong adverse effects with respect to the convergence during the solution stage.

A numerical solver for the electromagnetic magnetic voltage potential equation needs to address these problems that are generated by the large kernel of the matrix. To summarise, it should possess the following features:

• Iterative application

The size of the linear system of equations for a general induction heating application, where the global domain Ω is discretised, can attain several million unknown degrees of freedom. In addition, it is often not possible to apply direct solution techniques, due to the large demand for storage space. Furthermore, the large kernel of the curl-curl operator, in connection with the general bad conditioning of the linear operator can result in the failure of direct solvers. Therefore, it should be possible to apply the linear solver iteratively.

• Scalability Due to the above-mentioned large computational demand, it is natural that the linear solution stage ought to be done in parallel. The application time should scale inversely with the number of computational nodes, so that problems with increased computational demand can be solved efficiently.

• Cheap setup costs

The discrete bilinear form of the magnetic vector potential equation (II.80) depends on the electrical conductivity σ, which is temperature dependent, as well as the relative magnetic permeability µ r , which depends on the temperature and the applied magnetic field. In general, the weak coupling procedure results in a repeated application of the linear solver to a changing linear system of equations. A possible change of the computational domain, due to the introduction of movements, complicates this further. It follows, that the solver should be adaptable to changing material coefficients and cheap to setup, so that the impact of the changing material coefficients and discretisation has no large impact on the overall solution time.

• High reduction of numerical error

The induction heating model is based on generating a heating source term, based on the transient form of the magnetic vector potential. It follows that any error in the magnetic vector potential computation transports directly to the heating source computation. Any error is then redistributed to the computation of this vector potential, due to the coupling of the material quantities. Therefore, it is important that the numerical error of the solution can sufficiently be reduced. In addition, the computational cost should scale linearly with increased error reduction.

• Incorporation of the problem of the large null space

The solution technique should directly address the underlying numerical problems of the discretisation. Therefore, it ought to address the large kernel that is generated by the application of the curl-curl operator. The result should be an advantageous solution technique that generates conforming solutions.

In the following sections, an overview is given over some commonly applied classical solution techniques, as well as recently published multigrid algorithms that address the large kernel of the linear system of equations. Initially, the advantages and shortcomings of the classical solution techniques are detailed, including the context in which they have to be placed with regard to the subsequent multigrid techniques. Afterwards, a multigrid technique for the electromagnetic problem is presented. It is shown how it can be implemented in a parallel context, including a section on coarsening procedures, which is a large research topic for algebraic multigrid methods.

It is mentioned above that it should be avoided to use direct solution techniques, due to their large computational demand and problems with regard to non-definiteness. Alternative solution techniques that can be used are the family of iterative Krylov subspace solvers. Let r 0 be the initial residual of the computation and r n the n-th residual attained at the n-th application of the iterative solver. A classical Krylov subspace correction technique is the conjugated gradient (CG) method, proposed in [START_REF] Hestenes | Methods of Conjugate Gradients for Solving Linear Systems[END_REF] (cf. [START_REF] Weiss | A theoretical overview of Krylov subspace methods[END_REF]). Its repeated application creates a polynomial p n (A), defined by

r n = p n (A)r 0 , (III.2)
which minimises a quadratic function

f (x) = 1 2 x T Ax -b T x + c, (III.3)
which is based on the solution vector x, the linear system of equations A, the right hand side vector b and some constant c. Unfortunately, the convergence is not monotonic, which is why in general, the convergence shows an oscillating behaviour.

For exact arithmetic it is proven that the CG method will converge in at least N steps, where N is the dimension of the linear system of equations, if the matrix A is symmetric and positive definite. For non-symmetric matrices convergence can generally not be achieved. In this case, alternative methods need to be considered, like the bi-conjugated gradient method (BICG), presented in [START_REF] Fletcher | Conjugate Gradient Methods for Indefinite Systems[END_REF]. It introduces an additional polynomial pn (A T ) based on the transposed form of the linear operator that introduces a bi-orthogonality for the residual error. Due to the fact that it leads to large oscillations during the solution stage (see, e.g. [START_REF] Weiss | A theoretical overview of Krylov subspace methods[END_REF]), alternative techniques have been proposed, like the conjugate gradient squared (CGS) method in [START_REF] Sonneveld | CGS, A fast Lanczos-type solver for nonsymmetric linear systems[END_REF] or the stabilised BICG method (BICGStab) in [ [START_REF] Van Der Vorst | BI-CGSTAB: A fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems[END_REF]. The CGS method introduces a squared form of the initial matrix polynomial of the CG method, presented in equation (III.2). [START_REF] Kelley | Iterative Methods for Linear and Nonlinear Equations[END_REF] demonstrates that this increases the oscillation effects and can lead to instability and a loss of convergence. The BICGStab method introduces an additional polynomial Φ n (A), which is applied to equation (III.2). It is defined, so that it solves a one-dimensional optimisation problem in every iteration. The advantage is a generally smoother convergence, as shown by [START_REF] Weiss | A theoretical overview of Krylov subspace methods[END_REF] and [START_REF] Van Der Vorst | BI-CGSTAB: A fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems[END_REF], in comparison to the classical methods. Yet, its convergence is not monotonic. A comparable approach to the above-mentioned BICGStab method has been utilised in the generalised minimal residual (GMRES) method, which is proposed in [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]. It generates an n-dimensional orthogonal subspace that is increased in every iteration. In this method the residual r n is minimised for the full subspace in every iteration. This leads to smooth convergence with monotonic behaviour, i.e. the n+1-th residual is always smaller than or equal to the residual of the former iteration, so that ∥r n+1 ∥ ≤ ∥r n ∥.

(III.4)

The disadvantage of the GMRES method is the need to store an n-dimensional subspace that grows with every iteration. If convergence can be attained in a reasonable number of iterations it is more advantageous to employ the GMRES method over comparable solution techniques like the CG method.

III.2 Preconditioners

The iterative solution methods are, in general, not directly applied to the linear system of equations. Instead, a preconditioning operator is introduced that enhances the convergence of the solution procedure. Following [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF], the objective of a preconditioner is an increase of the robustness and efficiency of the solution method. It states that for most real world problems the quality of the preconditioning technique is more important than the choice of Krylov subspace solution method. The overall goal of the preconditioning technique is to approximate the inverse of the linear system of equations in a flexible, yet economical manner, so that the repeated application of the subspace correction technique leads to fast convergence. [START_REF] Kelley | Iterative Methods for Linear and Nonlinear Equations[END_REF] notes that this can be accomplished if the preconditioner clusters the eigenvalues of the linear system of equations close to one, which is the classical motivation for the introduction of a preconditioning operator. More recent developments in the field of electromagnetics show that this does not suffice [START_REF] Hiptmair | Auxiliary Space Preconditioning for Edge Elements[END_REF]. It shows that the correct replication of the vector components of the residual vector is essential, in order to steadily reduce the global error. The following sections introduce common preconditioners that are employed in the literature, showing the advantages and disadvantages with respect to an application in electromagnetics, serving as an introduction to the newer subspace splitting methods. These methods can be used in connection with multigrid methods, in order to effectively treat the large kernel of the linear system of equations, based on a correct reproduction of the underlying finite element spaces of the discrete approximations. In contrast to the above-mentioned classical preconditioners, this technique greatly reduces the effort of finding a numerical solution, so that the overall computational time is reduced, which opens up new possibilities for the numerical modelling of complex geometrical parts, featuring millions of degrees of freedom.

III.2.1 Classical preconditioners

The classical preconditioners that are used to approximate the inverse of the given linear system of equations are the Jacobian diagonal scaling preconditioner, the Gauss-Seidel method (GS), as well as the successive over-relaxation method (SOR) (see, e.g. [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]). Let D = diag(M) be the diagonal of the linear system of equations, so that the preconditioned matrix system can be rewritten using the Jacobian preconditioner P -1 as

P -1 J Mx = D -1 Mx = D -1 b. (III.5)
If the matrix M is diagonally dominant, the condition number of the transformed linear system of equations will be smaller than that of the original linear system of equations, so that convergence of Krylov methods can be achieved in a smaller number of iterations. The Gauss-Seidel preconditioner P GS , as well as the SOR precondtitioner P SOR are natural extensions of this method. These preconditioners are based on the matrix splitting M = L + D + R, where the matrix operators L and R denote the lower left triangular part and the upper right triangular part of matrix M. Both preconditioners are based on the approximative solution using the easily invertible matrix parts in either direct form (GS) or in weighted form using an overrellaxation parameter (SOR). Implementations can be found in either [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] or [START_REF] Kelley | Iterative Methods for Linear and Nonlinear Equations[END_REF]. For simple, well behaved problems, like the Laplacian or the the Poisson problem, convergence can be achieved comparatively fast, but the application does not assure fast convergence in any case. E.g [START_REF] Kelley | Iterative Methods for Linear and Nonlinear Equations[END_REF] notes that for general applications to numerical discretisations of partial differential equations "these preconditioners maybe somewhat useful, but should not be expected to have dramatic effects". [START_REF] Briggs | Introduction[END_REF] notes that the remaining residual error during an iterative solution is built from a low frequency error component e nl , as well as a high frequency components e nh . [START_REF] Briggs | Introduction[END_REF] shows that, in essence, the above-mentioned preconditioners reduce the high frequency parts, but have almost no effect on the low frequency parts. Therefore, even though convergence is assured at least at the N -th iteraton when exact arithmetic is used, in practice more elaborate techniques should be used.

Another particular family of classical preconditioners tries to approximate the inverse directly using multiplicative decompositioning techniques. Let L and R be triangular matrices, so that the invertible matrix M can be expressed by M = L Ū. The matrix inverse of M can now easily be expressed by the application of the inverse of each triangular matrix decomposition. Unfortunately, in practice, this LU-decomposition is computationally intensive and often results in dense decompositioning matrices [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. Therefore, an incomplete LU-decompositioning technique (ILU) is often employed that is based on a reduced fill-in, based on a given matrix fill-in ordering or a prescribed threshold tolerance. For symmetric positive definite matrices a Cholesky decomposition can be defined, so that the relationship L = U T is fulfilled, which results in storage savings. The classical application domain for this family of preconditioners are again the family of PDEs that are well-behaved and easily invertible. Yet, in the context of electromagnetism, the ill-conditioning of the linear system of equations results in adverse effects, during the solution stage. In practice, the application of the ILU-preconditioning technique can lead to destabilising behaviour, as presented in figure III.1b, which shows the convergence history of a preconditioned GMRES solver applied to the magnetic vector potential computation of the given crankshaft test case. The monotonicity of the GMRES method prevails, even though the application of an incomplete Cholesky type preconditioner (IC) leads to a severe flattening of the convergence rate. This can be attributed to the bad approximation qualities of the approximate inverse, which are enhanced by floating point inaccuracies. A simple Jacobian diagonal scaling preconditioner shows a better convergence rate, which is surpassed by a symmetric version of the SOR technique. Yet, for each of these classical preconditioning techniques it can be shown that the convergence rate flattens considerably, so that a repeated application during a solution phase is infeasible. The same figure details the results of an application of an auxiliary subspace correction technique that reduces the numerical error in connection with an algebraic multigrid method. It can be seen that convergence is achieved quickly, which can be attributed to the correct treatment of the large null space of the linear system of equations. The following sections detail the key concepts of this auxiliary subspace splitting method, including the multigrid method with its application to electromagnetic magnetic vector potential problems. The resulting auxiliary space multigrid method is easily applicable to standard edge element based discretisations and enables the efficient repeated solution of global electromagnetic problems, so that the modelling of large engineering simulations becomes feasible.

III.2.2 Stable decomposition of H 0 (curl, Ω)

The above-mentioned preconditioners are effective in reducing high frequency residual errors. The peculiarity of the curl-conforming finite element description, which conformingly discretises the large kernel, necessitates a more direct treatment of the particular error components. [START_REF] Hiptmair | Auxiliary space preconditioning in H0(curl; Omega)[END_REF]] introduces a space decompositioning technique that can be used to create effective preconditioners for every subspace of the approximating curl-conforming finite element space. Let H 0 (curl, Ω) be the square integrable curl-conforming inner product space, with vanishing surface values, defined through

H 0 (curl, Ω) := {u ∈ L 2 (Ω) : ∇ × u ∈ L 2 (Ω)}. (III.6)
The key idea of the stable vector space decomposition is that the vector space (III.6) can be decomposed into several corresponding Hilbert vector spaces. Let a generic inner product space be denoted by X i , so that x i ∈ X i , where i is chosen as counting index. The decomposition of any u ∈ H 0 (curl, Ω) into n orthogonal vector functions then follows as

u = x 1 + x 2 + • • • + x n , (III.7)
which can be rewritten in discrete form as

u h = x h1 + x h2 + • • • + x hn (III.8)
for the linear first order Nédélec finite elements u h ∈ E h ⊂ H 0 (curl, Ω) and the generic discrete Hilbert spaces x hi ∈ X hi . [START_REF] Hiptmair | Auxiliary Space Preconditioning for Edge Elements[END_REF] states that a preconditioner can be based on the discrete decomposition, if surjective bounded linear transfer operators Π hi can be generated for any space X hi , which define the map from the generic space to the vector finite element space

Π hi : X hi → -→ E h . (III.9)
Then, following from (III.8) a decomposition based preconditioner P -1 D can be generated as

P -1 D = n  i=1 Π hi X -1 hi Π T hi , (III.10)
where X hi is the norm inducing bilinear form on X hi and X -1 hi is an efficient preconditioner for the aforementioned bilinear form.

The choice of orthogonal subspace decomposition is quite arbitrary and depends on the nature of the underlying problem. For the problem of electromagnetics, an initial Helmholtz-decomposition has been proposed in ( [START_REF] Hiptmair | Multigrid for eddy current computation[END_REF], sec. 4). It is based on an L 2 -orthogonal decomposition of the curl-conforming space into a function space of the kernel of the curl operator, denoted by ker, and its complement. Let H 0 (Ω) be defined as the usual gradient space

H 0 (Ω) := {s ∈ L 2 (Ω) : ∇s ∈ L 2 (Ω)} (III.11)
and let (H 0 (Ω)) 3 be the corresponding space of vector functions s with respect to (III.11). The decomposition approach is based on the relationship (see [START_REF] Kolev | Parallel Auxiliary Sace AMG for H (curl) Problems. Technical report[END_REF], Thm. 3.1)

∀u ∈ H 0 (curl, Ω) ∃s ∈ (H 0 (Ω)) 3 : ∇ × s -∇ × u = 0. (III.12)
It is known that for the vector space H 0 (curl, Ω) the kernel of the differential curl operator can be represented by

H 0 (ker(curl), Ω) := {u ∈ H 0 (curl, Ω) : ∇ × u = 0} = ∇H 0 (Ω), (III.13)
so that following from (III.12) the two orthogonal subspaces can be identified as

H 0 (curl, Ω) = (H 0 (Ω)) 3 ⊕ ∇H 0 (Ω). (III.14)
Yet, [START_REF] Hiptmair | Nodal Auxiliary Space Preconiditioning in H(curl) and H(div) Spaces[END_REF] notes (cf. [START_REF] Hiptmair | Auxiliary Space Preconditioning for Edge Elements[END_REF]) that the orthogonal decomposition (III.14) is only valid in the context of a continuous formulation. [START_REF] Hiptmair | Auxiliary Space Preconditioning for Edge Elements[END_REF] mentions that a discrete decomposition is not valid, since local oscillatory functions can not be captured in the discrete decomposition. This contribution is identified as a high frequency contribution in [START_REF] Hiptmair | Nodal Auxiliary Space Preconiditioning in H(curl) and H(div) Spaces[END_REF], which will be denoted as ũh ∈ E h . This approach is consistent with the initial motivation of choosing Nédélec vector finite elements, as presented in chapter II. The main motivation for using the curl-conforming Nédélec edge vector elements is the correct reproduction of the underlying physical phenomena, which can not be captured by a nodal description.

Let N h and (N h ) 3 be the discrete representations of H 0 (Ω) and (H 0 (Ω)) 3 . The discrete version of s is identified in [START_REF] Hiptmair | Nodal Auxiliary Space Preconiditioning in H(curl) and H(div) Spaces[END_REF] in discretely projected form (see also [START_REF] Kolev | Parallel Auxiliary Sace AMG for H (curl) Problems. Technical report[END_REF] , Thm. 2.1) using the Nédélec interpolation operator N h as the sum over the n e edges e of the finite element discretisation, using the line integral over the edge tangent t e , multiplied by the corresponding local test function Ψ e as defined in equation (II.103) as

N h s h = ne  e=1 (  e s • t e de)Ψ e , ∀s h ∈ N 3 h . (III.15)
Similarly, the gradient of the scalar potential s can be defined in discrete form for all s h ∈ N h using the gradient operator G h . Finally, a stable decomposition of the form (III.8) can then be given using three discrete projection operators Π h1 = I, Π h2 = N h and Π h3 = G h as

u h = Π h1 x h1 + Π h2 x h2 + Π h3 x h3 = Iũ h + N h s h + G h s h .
The discrete decomposition based preconditioner (III.10) can then be defined as

P -1 D = Π h1 X -1 h1 Π T h1 + Π h2 X -1 h2 Π T h2 + Π h3 X -1 h3 Π T h3 (III.16) = X -1 h1 + N h X -1 h2 N T h + G h X -1 h3 G T h .
(III.17)

The bilinear form X h1 can be identified as the discretised form of the standard bilinear form of the magnetic vector potential problem (II.80), i.e. a(u h , v h ) for u h , v h ∈ E h . [START_REF] Hiptmair | Nodal Auxiliary Space Preconiditioning in H(curl) and H(div) Spaces[END_REF] notes that the error components are essentially due to high frequency components. Therefore, [START_REF] Hiptmair | Auxiliary Space Preconditioning for Edge Elements[END_REF] states that for X h1 a simple smoothing technique, like introduced in subsection (III.2.1), can be chosen as preconditioning method X -1 h1 . The bilinear forms X h2 and X h3 can be chosen more freely. The definition of these bilinear forms is essentially based on the choice of auxiliary discretisation method. E.g. the bilinear forms can be defined on an auxiliary regular mesh, so that the auxiliary problems are amenable for a treatment with geometric multigrid methods. Alternatively, a discretisation can be based on the existing nodal description of the finite element discretisation, with the help of suitable projection operators. A detailed overview over common methods is given in [START_REF] Kolev | Parallel Auxiliary Sace AMG for H (curl) Problems. Technical report[END_REF]. In the following, some common methods are detailed in order to illustrate the versatility of the decomposition based preconditioning technique.

III.2.2.1 Auxiliary mesh method

The auxiliary mesh method has been proposed in [START_REF] Hiptmair | Auxiliary space preconditioning in H0(curl; Omega)[END_REF]. The initial problem is defined on a locally refined and unstructured finite element discretisation T , discretised using linear edge elements E h . [START_REF] Hiptmair | Auxiliary space preconditioning in H0(curl; Omega)[END_REF]] introduces a regular auxiliary mesh T a , discretised using vector type finite elements of second order Ẽh [START_REF] Nedelec | A new Family of Mixed Finite Elements in R3[END_REF], with two degrees of freedom per edge. In this context, the interpolation operator Π h2 is identified as interpolation operator from a regular auxiliary mesh to the initial mesh

Π h2 = T : Ẽh → -→ E h .
(III.18)

The discrete preconditioner (III.10) is then defined as

P -1 D =X -1 h1 +T h X -1 h2 T T h +G h X -1 h3 G T h , (III.19)
where X h1 describes the bilinear form on T using E h and X h2 is the bilinear form on T a using Ẽh . X h3 is chosen as the variationally equivalent Poisson restriction of the original bilinear form

X h3 = G h X h1 G T h .
(III.20) [START_REF] Hiptmair | Auxiliary space preconditioning in H0(curl; Omega)[END_REF] proposes a Gauss-Seidel smoothing technique for the preconditioners X -1 h1 and X -1 h3 , whereas X -1 h2 is defined using a geometric multigrid method on the auxiliary mesh T a . Details for multigrid methods applied to the curl-conforming finite element vector space can be found in [START_REF] Hiptmair | Multigrid Method for Maxwell's Equations[END_REF] and [START_REF] Arnold | Multigrid in H(div) and H(curl)[END_REF]. For the auxiliary mesh method, [START_REF] Hiptmair | Auxiliary space preconditioning in H0(curl; Omega)[END_REF] gives some results regarding two-dimensional sample test problems, showing good convergence results when the preconditioner is used in connection with a CG solver. Furthermore, [START_REF] Hiptmair | Auxiliary space preconditioning in H0(curl; Omega)[END_REF] states that the considerations regarding the auxiliary finite element space can be relaxed, so that in a practical application the auxiliary space can be chosen as E h , but concluding that problems with strongly varying material coefficients might need entirely different bilinear forms on the auxiliary space. The auxiliary mesh preconditioning approach has also been studied by [START_REF] Kolev | H (curl) auxiliary mesh preconditioning[END_REF], which gives results for several numerical experiments. [START_REF] Kolev | H (curl) auxiliary mesh preconditioning[END_REF] shows that the approach is, in general, very efficient, but depends on the chosen auxiliary mesh, stating that, in comparison to matching auxiliary meshes, non-matching auxiliary meshes lead to a reduction in the convergence rate. In addition, it is mentioned that, in comparison to an additive application, a multiplicative application of the preconditioners leads to an increase in the convergence rate.

III.2.2.2 Auxiliary nodal space method with discrete elliptic operators

The auxiliary nodal space method is based on the assumption that the auxiliary space can directly be discretised using the existing nodes of the finite element mesh. The decomposition based preconditioner can then directly be chosen as (III.17), with the given transfer operators Π h1 = I, Π h2 = N h and Π h3 = G h . Again, operator X h1 can be identified as the bilinear form on the initial edge finite element mesh, with suitable smoothing preconditioner X -1 h1 . In the discrete elliptic approach, both bilinear forms X h2 and X h3 are explicitly discretised using a weak formulation. [START_REF] Hiptmair | Auxiliary Space Preconditioning for Edge Elements[END_REF] proposes the forms

X h2 =  Ω α∇u h ∇v h dΩ, ∀u h , v h ∈ (N h ) 3 (III.21) X h3 =  Ω τ ∇u h ∇v h dΩ, ∀u h , v h ∈ N h , (III.22)
whereas [Kolev and Vassilevski, 2006a] proposes to use the original bilinear form, so that

X h2 =  Ω τ u h v h dΩ+  Ω α∇u h ∇v h dΩ, ∀u h , v h ∈ (N h ) 3 (III.23) X h3 =  Ω τ ∇u h ∇v h dΩ, ∀u h , v h ∈ N h . (III.24)
The parameter α and τ depend on the magnetic permeability µ = µ 0 µ r , the electrical conductivity σ and the time step size ∆t. In this work, they can be defined as After conducting numerical experiments, [START_REF] Hiptmair | Auxiliary Space Preconditioning for Edge Elements[END_REF] concludes that the preconditioner is robust with a convergence that is mesh independent. The use of this auxiliary space preconditioner, in connection with a one level algebraic multigrid solver for the elliptic weak forms, leads to a rapid convergence in connection with a CG solver. [START_REF] Kolev | Parallel Auxiliary Sace AMG for H (curl) Problems. Technical report[END_REF] (cf. [Kolev and Vassilevski, 2006a]) gives additional results regarding large coefficient jumps, including a semi definite problem with vanishing electrical conductivity, in a highly parallel setting on up to 1024 processors for problems up to approximately 70-100 million degrees of freedom. It is shown that the method is efficient and scalable in connection with a one level algebraic multigrid preconditioner for both bilinear forms. In addition, it is stated that a variationally equivalent formulation of (III.23) and (III.24) shows an even better performance. It is furthermore shown that the auxiliary space preconditioner behaves qualitatively like a comparable Laplacian example problem that is solved using standard nodal algebraic multigrid. [START_REF] Kolev | Parallel Auxiliary Sace AMG for H (curl) Problems. Technical report[END_REF] conjectures that any advancement in standard algebraic multigrid method should be applicable to the auxiliary space preconditioning technique.

III.2.2.3 Auxiliary nodal space method with variationally equivalent elliptic operators

The variationally equivalent auxiliary nodal space method has been proposed in [Kolev and Vassilevski, 2006a]. Here, the auxiliary operators of the elliptic forms in equation (III.17) are specified using variationally equivalent bilinear forms. The transfer operators remain as Π h1 = I, Π h2 = N h and Π h3 = G h . The bilinear form X h1 is the bilinear form on the initial edge finite element mesh, so that a smoothing preconditioner X -1 h1 can be utilised. [Kolev and Vassilevski, 2006a] notes that the elliptic operators can be generated by a weak Galerkin projection approach, using the existing transfer operators N h and G h , so that X h2

X h2 = N T h X h1 N h , (III.27)
whereas X h3 can simply be defined as

X h3 = G T h X h1 G h . (III.28)
Again, the preconditioners for each of the resulting bilinear form can, in principle, be chosen arbitrarily. [Kolev and Vassilevski, 2006a] notes that the forms provided in equations (III.27) and (III.28) are essentially equivalent to the formulation provided in equations (III.23) and (III.24), so that one-cycle multigrid methods are sufficient to approximate the inverse forms. Its application is simple, since the weak formulations can be generated automatically, from the existing bilinear form of the original problem and both linear transfer operators N h and G h that have to be created for the transformation. [Kolev and Vassilevski, 2006a] recommends formulation (III.27) and (III.28), since the discretely formed elliptic operators, defined in equations (III.23) and (III.24), are not variationally equivalent to the initial bilinear form, so that convergence is not guaranteed.

III.2.2.4 Discretisation of N h and G h

The discrete gradient matrix maps from vertex space N h to the curl-conforming finite element space E h , so that

G h : N h → -→ E h . (III.29)
For the space of linear Nédélec vector finite elements, the discrete version is simple to generate. Following [START_REF] Hiptmair | Auxiliary Space Preconditioning for Edge Elements[END_REF], in this configuration the gradient can be approximated for each edge e as the vertex difference of each edge endpoint vertex i and j. The dimensions of G h then follows from the number of edges and vertices in the finite element mesh. The matrix G h is very sparse with only two non-zero entries per row. Let G h [i, j] denote the entry for the i-th row and j-th column of the discrete operator G h . Its row-wise definition then follows for each corresponding edge e as

G h [e, i] = -1 (III.30) G h [e, j] = 1. (III.31)
The number of rows equals the number of edge degrees of freedom n e , whereas the number of columns equals the number of vertex nodes n v , so that its structure can be visualised as

G h =       -1 1 -1 1 . . . . . . -1 1       (III.32) n v n e .
According to [START_REF] Hiptmair | Auxiliary Space Preconditioning for Edge Elements[END_REF], the transfer operator N h can be defined in a similar manner, by noting that it details an equally weighted transfer of nodal values to the edge. It has to be accomplished for three components in R 3 , so that it can be directly based on the gradient operator (III.29), which defines an unweighted vertex difference. Its discrete form can be given in block form as

N h = 1 2  N x h N y h N z h  (III.33) 3 × n v n e ,
where 1 2 is the weighting factor for the two edge vertices and N x h , N y h and N z h are block matrices with the same sparsity structure as G h (see [START_REF] Hiptmair | Nodal Auxiliary Space Preconiditioning in H(curl) and H(div) Spaces[END_REF], eq. 7.4, [START_REF] Kolev | Some experience with a H1-based auxiliary space AMG for H (curl) problems[END_REF], sec. 5 and [START_REF] Kolev | Parallel Auxiliary Sace AMG for H (curl) Problems. Technical report[END_REF], sec. 5). The number of rows equals the number of edge degrees of freedom n e , whereas the number of columns equals the number of vertices n v multiplied by the dimension of the vector space, i.e. three for R 3 . The block matrices can be formed after application of the gradient matrix G h to the vertex coordinates x, y and z in R 3 . The i-th and j-th column entry for row e equals the corresponding row of the matrix vector product of the gradient operator and the vertex nodes. The first block matrix N x h can be defined as

N x h [e, i] =(G h x)[e] (III.34) N x h [e, j] =(G h x)[e],
(III.35) whereas N x h and N x h are defined as

N y h [e, i] =(G h y)[e] (III.36) N y h [e, j] =(G h y)[e] (III.37)
and

N z h [e, i] =(G h z)[e] (III.38) N z h [e, j] =(G h z)[e].
(III.39)

It can be noted that this adjacency information is readily available for edge finite element applications. Therefore, the implementation costs are very low.

III.3 Multigrid methods

Let n denote the current iteration of the numerical solution method. [START_REF] Briggs | Introduction[END_REF] states that the residual error e n for the n-th iteration can be split into two components as e n = e nh + e nl .

(III.40)

The first component e nh is denoted as high frequency component, whereas the second component e nl is described as low frequency component. The denomination stems from the observation that a smoothing operator S, like the standard Jacobian or Gauss-Seidel method (see subsection III.2.1) can very effectively reduce a given component e nh , but fail to reduce e nl [START_REF] Brandt | Guide to multigrid development[END_REF]. The behaviour is mesh dependent, so that an increase or decrease of the characteristic element size changes the behaviour of the smoothing operator with respect to the numerical error. If the residual is restricted to a coarser discretisation, the error component e nl starts to behave more like a high frequency component e nh , so that smoothing becomes more effective.

A multigrid method uses this effect to rapidly decrease the residual error. The most basic form is a V-cycle [START_REF] Barth | Multiscale and multiresolution methods: Theory and applications[END_REF]. The number of coarsening levels can be adapted to the underlying problem and depends on the available computing resources. In the following, the key ideas of the multigrid method are demonstrated using a three-step scheme including two fine discretisations. Figure III.2a gives an overview of this three-level multigrid discretisation. The application ordering of the operators of the corresponding standard V-cycle is depicted in III.2b. Initially, a smoothing operator S 1 is applied to an initial discretisation, which is denoted in discrete matrix form as A 1 , with index 1 denoting the first level. This results in a decrease of the high-level error component. Afterwards, a restriction operator R 12 is used to restrict the computed residual to a coarser level, which can be denoted as A 2 . The indices denote the direction of the transformation from fine grid 1 to coarse grid 2. As stated above, the remaining low frequency error of the residual component becomes a high-level component on the coarse grid, so that a second smoothing operation S 2 can be applied. This recursive procedure can be done until the coarsest level A m has been reached at the m-th level. Afterwards, the restricted residual has to be interpolated up to the fine initial discretisation. This is achieved by using an interpolation operator P jk that prolongates the results of the j-th level to a finer level k.

In the classical multigrid scheme, the discrete operators A 1 , . . . , A m and the respective restriction and prolongation operators are defined using nested discretisations. Restrictions and prolongations are particularly easy to accomplish, if the discretisation is based on a regular mesh. Non-regular finite element meshes necessitate a treatment using special meshing software, in order to create nested finite element sub-meshes. Unfortunately, it is particularly difficult to create these for higher dimensional problems and finite element discretisations involving complex geometries. Example implementations for electromagnetic problems are presented in [START_REF] Hiptmair | Multigrid Method for Maxwell's Equations[END_REF], [START_REF] Hiptmair | Multigrid for eddy current computation[END_REF] and [START_REF] Mifune | A Fast Solver for FEM Analyses Using the Parallelized Algebraic Multigrid Method[END_REF]. The method presented in [START_REF] Hiptmair | Multigrid Method for Maxwell's Equations[END_REF] and [START_REF] Hiptmair | Multigrid for eddy current computation[END_REF] is based on a Helmholtz-type orthogonal subspace decomposition, i.e. it is based on the decomposition of the vector field F into an irrotational and a rotational part, so that

F = F irr  ∇×F irr =0 + F rot  ∇•Frot=0 , (III.41)
which results in the continuous decomposition of the kernel of the differential curloperator and its complement, demonstrated in the equation (III.14). The multigrid schemes are based on nodal multigrid techniques that conform to the initial edge finite element description. These methods form the basis of the auxiliary subspace correction schemes. [START_REF] Hiptmair | Multigrid Method for Maxwell's Equations[END_REF] demonstrates the good convergence behaviour for this multigrid scheme for problems with homogeneous materials. It is furthermore demonstrated that it is not heavily influenced by the time step size of an implicit discretisation. Yet, [START_REF] Hiptmair | Multigrid Method for Maxwell's Equations[END_REF] reports that despite some success for problems with non-homogeneous materials, the convergence rate deteriorates if the number of non-homogeneous domains is increased. An enhancement to the classical multigrid method is the algebraic multigrid method [START_REF] Stüben | A review of algebraic multigrid[END_REF]. In contrast to the geometrical multigrid method, all coarsening operators are based on a coarsening of the initial matrix representation A 1 . The resulting discretisations are comparable to the two-dimensional example discretisation, visible in figure III.2a. The initial mesh is non-regular. The coarser meshes are defined by removing certain edges and nodes. A possible coarsening method that is based on the matrix discretisation A 1 uses the same approach to generate coarser discretisation levels. This means that discretisations like A 2 or A 3 in figure III.2a can directly be specified using recursive relationships with respect to A 1 . The application is, therefore, not dependent on an outside coarsening of the underlying discretisation. Therefore, the resulting method is very flexible and adaptable, since there is no need for special meshing software. Example applications for electromagnetic problems are presented in [START_REF] Mifune | A Fast Solver for FEM Analyses Using the Parallelized Algebraic Multigrid Method[END_REF], [START_REF] Lee | A Novel Algebraic Approach for Maxwell's Equations[END_REF] and [START_REF] Kolev | Parallel Auxiliary Sace AMG for H (curl) Problems. Technical report[END_REF].

The method demonstrated in [START_REF] Mifune | A Fast Solver for FEM Analyses Using the Parallelized Algebraic Multigrid Method[END_REF] is a classical nodal multigrid scheme for a nodal finite element description, showing good parallel efficiency and generally good convergence for a simple two-dimensional example domain in structured and unstructured form. But it has to be noted that the discretisation problems of the nodal finite element method hinder the application of this method in more complex cases.

[ [START_REF] Kolev | Parallel Auxiliary Sace AMG for H (curl) Problems. Technical report[END_REF] demonstrates the application of the auxiliary subspace decomposition technique using the efficient algebraic multigrid implementation presented in [START_REF] Henson | BoomerAMG : A parallel algebraic multigrid solver and preconditioner[END_REF], showing an excellent scalability and numerical efficiency for large scale applications. [START_REF] Kolev | Parallel Auxiliary Sace AMG for H (curl) Problems. Technical report[END_REF] demonstrates that it suffices to use a one-level V-cycle type multigrid method to achieve a high convergence rate.

[ [START_REF] Lee | A Novel Algebraic Approach for Maxwell's Equations[END_REF] proposes a Helmholtz-type decomposition based multigrid scheme using a tree and co-tree based decomposition of an edge finite element discretisation. Several test cases are presented, including examples with varying coefficients. The viability of the method is demonstrated, but the splitting into tree and co-tree spaces results in a coupled linear system of equations in block form, resulting in an implementation that is not as simple as the auxiliary subspace correction technique used in [START_REF] Kolev | Parallel Auxiliary Sace AMG for H (curl) Problems. Technical report[END_REF].

It has to be noted that the coarsening procedure depends on the underlying problem and the available computing resources. In the following an overview of some available coarsening schemes is given. The choice of coarsening has an impact on the convergence behaviour of the multigrid method, the setup costs, but also on the parallelisability.

The multigrid method is especially efficient for solving problems of elliptic type. Therefore, these methods are particularly efficient for solving the preconditioning technique presented in subsection III.2.2, which is based on a stable subspace decomposition using elliptic forms. The application of the algebraic multigrid method results in very versatile and flexible preconditioners that can increase the convergence rate of standard Krylov subspace methods.

III.3.1 Coarsening techniques

An application of the algebraic multigrid method necessitates a coarsening of the fine matrix operators. In the following, some common coarsening schemes are detailed, which can be used in connection with the auxiliary subspace correction technique. Each coarsening algorithm is briefly explained and the key ideas that lead to its development are presented. Afterwards, the advantages and shortcomings of each method are highlighted.

III.3.1.1 Requirements

Coarsening algorithms are interpolation and reduction algorithms relating some nodes F of a fine mesh to a coarse nodal subset C, so that the reduction relationship to the coarser level describes the relationship

G N = C ∪ F, (III.42)
where the full set of nodes in the graph of the initial level is denoted as G N . The nodal subset F builds the set of interpolatable nodes, which are the complementary part for the coarse set C on the coarser level k, so that

C ∩ F = ∅. (III.43)
Following [START_REF] Henson | BoomerAMG : A parallel algebraic multigrid solver and preconditioner[END_REF], the interpolation describes the numerical error e i for each fine point i ∈ F as a function of coarse grid points, so that the error can be reproduced on the finer grid using the coarse interpolation nodes of the subset C i ⊂ C as

e i =  j∈C i w ij • e j , (III.44)
where w ij is a weighting factor that describes the influence of the coarse point with index i on the fine point with index j and e j is the numerical error at the j-th coarse point. A stabilised interpolation formula for the weightings w ij is presented in [START_REF] Henson | BoomerAMG : A parallel algebraic multigrid solver and preconditioner[END_REF]. [START_REF] Henson | BoomerAMG : A parallel algebraic multigrid solver and preconditioner[END_REF] states that the modified interpolation weights can be defined as

w ij = - 1 a ii +  k∈(D w i ∪F i ) a ik   a ij +  k∈D s i \F i a ik âkj  m∈C i âkm   . (III.45)
Here, the values a ii , a ik , a ij and a ik are entries of the discretisation matrix A on the corresponding level and âkj and âkm are defined as The nodal subsets D w i describes the points connecting to i that influence it only weakly, whereas D s i describes the points that strongly influence i. Both of these subsets do not coincide with the coarse points C, i.e. D w i ̸ ⊂ C and D s i ̸ ⊂ C. This modified formula for generating weighting values slightly deviates from the classical interpolation formula and has been derived by [START_REF] Henson | BoomerAMG : A parallel algebraic multigrid solver and preconditioner[END_REF] to avoid extremely large interpolation weights, which are sometimes generated in the classical formula and which might lead to numerical errors and divergence.

âkj =  0, if sgn(a kj ) =
The defining operation for each type of algebraic multigrid method is the choice of coarsening. For each coarsening scheme a heuristic is described, which forms the basis of the algorithm that divides G N into F and C.

According to [START_REF] Sterck | Reducing Complexity in Parallel Algebraic Multigrid Preconditioners[END_REF] the heuristics of common coarsening schemes can be defined as H1: If j strongly influences i ∈ F then either j ∈ C or j strongly (III.48) depends on a k ∈ C that also strongly influences i. H2: C should form a maximal independent set [START_REF] Tutte | Graph Theory[END_REF] in G N . (III.49) It is mentioned that H1 ensures adequate interpolation, so that the convergence rate of the algebraic multigrid method is constant and bounded away from one, whereas H2 is chosen, so that there is a sufficient number of coarse points to accurately interpolate, but not more than necessary. In that regard the term strongly denotes a connection for which

-a ij ≥ α t max(-a ik ), ∀k ̸ = i (III.50)
is valid. The variable α t is a supplied threshold value that reduces the storage demand of the discrete operators. Its value is chosen between 0 < α t ≤ 1. A large number of connections generally leads to an increase in the convergence rate, but at the same time, to an increase in the storage and application cost. Additionally, two ratios can be defined, which describe the behaviour of an algebraic multigrid method. The operator complexity c is defined as the quotient of the sum of non-zero values, defined as nnz, of all discretisations on all m discretisation levels, with respect to the initial level. The application of the multigrid method in the context of this work is related to a one-cycle multigrid method, so that the operator complexity reduces to

c = nnz(A 1 ) + nnz(A 2 ) nnz(A 1 ) . (III.51)
The stencil s is a coefficient value expressing the average number of non-zero values per row for each respective discretisation level. It is motivated by noting that a large stencil size leads to an increase in setup time and application cost, since the discretisation becomes more dense. [Meier Yang, 2006] notes that this is a common problem for multiple coarsening levels, so that the coarsest level possesses a large stencil size. In the context of the auxiliary subspace preconditioning technique, this value can be neglected, since only one coarse level is formed. Following is a brief overview of common coarsening methods that have been used in this work to create coarse discretisations. A more general overview can be found in [START_REF] Henson | BoomerAMG : A parallel algebraic multigrid solver and preconditioner[END_REF] and [START_REF] Sterck | Reducing Complexity in Parallel Algebraic Multigrid Preconditioners[END_REF].

III.3.1.2 RS, RS2 and RS3

The Ruge-Stüben (RS) method has been proposed in [START_REF] Ruge | Algebraic multigrid[END_REF]. This method is based on a variable λ i that counts the number of points j, which are strongly influenced by i. From all the nodes of the nodal set with the maximum connections, one node is randomly picked and designated as coarse node in C. All nodes that depend on i are included in the set F . These interpolated points are connected to neighbour points, denoted as k. If these neighbour points strongly influence the points in the set F , each corresponding value λ k is incremented. Afterwards, one node will be designated as additional coarse point and subsequently added to the set of coarse points C. The algorithm is then recursively applied, until every node has been assigned. [START_REF] Sterck | Reducing Complexity in Parallel Algebraic Multigrid Preconditioners[END_REF] states that counterexamples can be generated for which points i, j ∈ F can be found that invalidate the condition H1. For these cases correct interpolation is, therefore, not assured. A remedy is the second stage Ruge-Stüben (RS2) method that adds a second pass over all nodes in F . For every node i ∈ F it is checked, whether there exists a strong dependency between i and a different node j ∈ F that is not related to a common coarse point k ∈ C. If such a pair is found, one node will be assigned to the set C, so that a correct fine grid interpolation is assured.

The RS and RS2 methods are inherently sequential. A parallelisation can be implemented by dividing the nodal subset G N into smaller subsets, followed by a local application of RS2. This results in inconsistencies for the nodes on the processor boundaries. As a remedy, [START_REF] Henson | BoomerAMG : A parallel algebraic multigrid solver and preconditioner[END_REF] proposes a third pass Ruge-Stüben (RS3) method by simply applying the second stage of the RS2 algorithm to every boundary node, so that correct C sets are produced.

The RS method enforces neither H1 nor H2, whereas RS2 and RS3 enforce H1, but not H2. But it can be noted that in all cases H2 is approximatively fulfilled.

III.3.1.3 CLJP

The Cleary-Luby-Jones-Plassman (CLJP) algorithm is proposed in [START_REF] Cleary | Coarse-Grid Selection for Parallel Algebraic Multigrid[END_REF]. Let rand(0,1) be a random number in the range [0, 1]. The CLJP algorithm initially bases the coarse grid selection on a randomised definition of the number of influencing points, in order to ensure uniqueness of the set of maximal connecting nodes, replacing λ i by λi = λ i + rand(0, 1).

(III.52)

As noted by [START_REF] Sterck | Reducing Complexity in Parallel Algebraic Multigrid Preconditioners[END_REF], the CLJP-algorithm possesses local maximal measures, which contrasts to the global maximal measures of the original RS method. It follows that the coarsening points can be selected in parallel, which increases the scalability with respect to the sequential procedure of the RS method. [START_REF] Sterck | Reducing Complexity in Parallel Algebraic Multigrid Preconditioners[END_REF] notes that the advantage of the CLJP method is that it can be efficiently parallelised, while retaining the original scalability of the RS coarsening algorithm, but it is also mentioned that due to the randomness of the selection algorithm, fine grid information can be lost by the coarsening procedure, which might lead to a degradation of convergence, even though the algorithm might still be more advantageous for unstructured problems than the classical RS methods. The advantage of the CLJP method is that the coarsening is only dependent on the initial assignment of weighting factors. It follows that the coarsening structure is not influenced by the number of processors and that any communication is restricted to the communication between boundary sets.

III.3.1.4 Falgout

The Falgout scheme is a simple extension of the RS and CLJP algorithms and has been presented in [START_REF] Henson | BoomerAMG : A parallel algebraic multigrid solver and preconditioner[END_REF]. Initially, the parallel application of a node distribution algorithm results in the assignment of the untreated nodes to each local processor. For each of these nodes an RS2 method is performed to generate a locally conforming coarse set. Since the boundary nodes might still possess dependencies that are not treated by a correct coarsening, an initial first independent set D of coarse points C is defined locally. Every coarse point i in C is included in D, if it is entirely local for each processor and not adjacent to any boundary point. From this set of starting values, a second CLJP pass is performed to generate coarse points on the processor boundaries. [START_REF] Henson | BoomerAMG : A parallel algebraic multigrid solver and preconditioner[END_REF] notes that the quality of the interior coarsening is generally very good for RS methods, whereas it is comparatively worse for the CLJP method. The application of the Falgout scheme leads therefore to a good quality of the coarsening for the interior domains, while ensuring a correct treatment of the boundary nodes. [START_REF] Henson | BoomerAMG : A parallel algebraic multigrid solver and preconditioner[END_REF] states that the Falgout scheme can be understood as a CLJP scheme with a different initial independent set D. The disadvantage is stated as the large growth in complexities for three-dimensional problems, which is a common problem for the Falgout, CLJP, as well as for the RS method. It is furthermore noted that, if the grid structure is regular, the Falgout algorithm is in general more efficient with much lower complexities than the RS and CLJP methods. This advantage can not be shown for irregular grids. For a problem on an irregular grid it is demonstrated that the complexities of both the CLJP and the Falgout method are comparable so that the solution times are almost equal. For a three-dimensional Laplacian problem with varying material coefficients it is demonstrated that the Falgout method is slightly faster than the CLJP method, with complexities that are approximately three times smaller.

III.3.1.5 PMIS

The parallel modified independent set (PMIS) method is proposed in [START_REF] Sterck | Reducing Complexity in Parallel Algebraic Multigrid Preconditioners[END_REF], which notes that it can be understood as a simplified version of the CLJP algorithm. The goal is to create a fully parallel algorithm that minimises message passing in that it only necessitates an interaction between neighbouring processors. It is furthermore independent on the number of processors, leading to the same coarsening for different parallel configurations. The heuristic H1 is ignored and replaced by H1 ′ , which is identified as H1 ′ : There exists at least one point j ∈ C with a (III.53) strong connection to i for every point i ∈ F .

The technical details of the implementation of the PMIS algorithm are presented in [START_REF] Sterck | Reducing Complexity in Parallel Algebraic Multigrid Preconditioners[END_REF]. Three different categories of nodal set are defined, which can be identified as the coarse nodes C, the fine nodes F and the unassigned nodes R. The iterative algorithm proceeds until all remaining nodes of the nodal set R have been assigned to either F or C. In the beginning R = G n and the relationship (III.50) is used to identify the initial set of fine points F , by noting that the nodes that do not strongly influence any other nodes must belong to F . The first independent set that is assigned to C is chosen as a random initial set, weighted in favour of strongly influencing points, using relationship (III.52), i.e. for the remaining nodes R every node i is assigned to the set C for which the relationship λi ≥ λj , ∀j connecting to i (III.54) is fulfilled. Afterwards, every point j that is found to remain in R, which is strongly influenced by a new point i ∈ C, is included into F . The algorithm then proceeds for the remaining points R up until every node has been assigned, either to C or F . It follows that the final independent set C is a maximal independent set, so that condition H2, defined in (III.49), is enforced. [START_REF] Sterck | Reducing Complexity in Parallel Algebraic Multigrid Preconditioners[END_REF] notes that for example test cases, PMIS can be two to three times faster than both the CLJP and the Falgout schemes with operator complexities of less than half the size. But it is also noted that the aggressive coarsening leads to a reduction in the convergence rate.

III.3.1.6 HMIS

The hybrid modified independent set (HMIS) method is proposed in [START_REF] Sterck | Reducing Complexity in Parallel Algebraic Multigrid Preconditioners[END_REF]. Like the Falgout algorithm it is a hybrid method. It is based on an initial local treatment using the RS scheme, i.e. on all processors the distributed nodes are divided into local sets F and C using the classical RS algorithm. The nodes i ∈ C that are not on an internal boundary are chosen as first initial set C. The remaining boundary nodes form the set of remaining nodes R. These boundary nodes are then treated in a following step by using the PMIS method, which is used instead of a CLJP stage for the Falgout scheme. The strength of the HMIS method is a reduction of the operator complexity and a decrease in setup costs. The disadvantage is a reduction of the convergence rate. It must be noted that it enforces condition H1 ′ , whereas condition H2 is not strictly enforced due to the application of the RS method for the interior nodes.

[ [START_REF] Sterck | Reducing Complexity in Parallel Algebraic Multigrid Preconditioners[END_REF] notes that for a three-dimensional Laplacian problem on a regular domain, both the PMIS, as well as the HMIS method have much smaller complexities and reduced stencil sizes than the CLJP and Falgout coarsening methods, with the complexity being around six times smaller than for the CLJP scheme. In addition, it is demonstrated that the PMIS and HMIS algorithms lead to algebraic multigrid operators necessitating more iterations to converge, but it is also shown that the reduction of the complexity leads to an overall reduction of computational time. In addition, it is shown that the convergence rate of the PMIS and HMIS method degrades, if the number of processors is increased, which is not a large problem for the Falgout and CLJP method. For an irregular description of the Laplacian problem it is shown that the complexities for the PMIS and HMIS algorithms are around three times lower than for corresponding Falgout and CLJP schemes. Again, an increase in the number of processors leads to a decrease in the convergence rate, but the net effect is a reduction of overall solution time. After introducing spatially varying material coefficients into the above-mentioned Laplacian formulation it is noted that the convergence rate severely degrades. The number of iterations for a sequential application for both the PMIS and HMIS methods is around seven times larger than for corresponding CLJP and Falgout schemes. It is shown that for a highly parallel test case with 1000 processors the convergence rate deteriorates further, so that around 40 times more iterations are necessary to achieve convergence.

III.3.1.7 ECGC

An important aspect in parallel coarsening is the correct treatment of boundary nodes. [Griebel et al., 2006a] states that the classical schemes necessitate message passing between distinct boundaries, to ensure a correct interpolation by coarse grid points. The addition of coarse points on the processor boundaries, like for the RS3 scheme results in an increase of operator complexity. For the CLJP algorithm this even leads to an increase of coarse points in the inside of the domain. Therefore, [Griebel et al., 2006a] proposes the coarse grid correction (CGC) method, which allows a parallel coarsening without the need of direct communication during the initial coarsening phase. The key idea follows from the observation that the classical RS scheme leads to distinct coarsening results, if different starting points are chosen. This allows to create several initial distinct coarsenings on each local processor domain. Afterwards, a weighted graph of these multiply defined coarsenings can be created that can be used to identify the local subsets that best approximate the global coarsening, i.e. the local coarsenings are picked that best match the neighbour coarsenings. If any remaining points i ∈ F exist that are strongly coupled to a node j ∈ F , they are associated to C. The resulting algorithm has a slightly higher setup cost than the alternative methods, since many initial coarse sets need to be generated. But it is mentioned by [Griebel et al., 2006a] that the cost of the additional local coarsenings is negligible, since the decrease of the operator complexity leads to time savings in the solution stage.

The drawback of the CGC algorithm is that it can not handle arbitrary small coarsenings, i.e. the algorithm fails, if the number of nodes is relatively small compared to the number of processors. [START_REF] Griebel | Coarse grid classification-Part II: Automatic coarse grid agglomeration for parallel AMG[END_REF] mentions that this leads to smallest coarsenings that scale with the number of cpus n proc , so that the coarsest node size N behaves like

N = O(n proc ).
(III.55)

The extended coarse grid classification (ECGC) algorithm remedies this problem. It is proposed in [START_REF] Griebel | Coarse grid classification-Part II: Automatic coarse grid agglomeration for parallel AMG[END_REF]. The ECGC method extends the original CGC method, by replacing local coarsenings that can be interpolated by outside grid points using empty grids in the global selection step. This allows arbitrary levels of coarsenings on the global level up to a single point. [START_REF] Griebel | Coarse grid classification-Part II: Automatic coarse grid agglomeration for parallel AMG[END_REF] provides several example results for the CGC and ECGC method in comparison to a Falgout scheme. It shows that the convergence factors and solution times of the CGC and the ECGC method are comparable to the Falgout method. It is furthermore demonstrated that the ECGC method leads to a reduction in setup times for highly parallel test cases, with respect to the CGC method. It is shown that the setup times are comparable to the setup times of the Falgout scheme, but it is mentioned that the ECGC method leads to the fastest setup times for highly parallel test cases. The operator complexities for the ECGC method for highly parallel test cases are slightly smaller than the operator complexities for comparable configurations involving the Falgout method, which leads to a reduction of solution time.

III.4 Application

The following section demonstrates the efficient solution of the magnetic vector potential problem for the example test case of an automotive crankshaft, embedded in an enclosing air domain. The source current is applied using a closed inductor of complex geometry. The finite element mesh includes around 3 million finite elements, resulting in around 3.5 million edge degrees of freedom. The sparse linear system of equations consists of approximately 55 million non-zero values. The impact on the different operator settings are detailed for common operator choices for the algebraic multigrid solver, which addresses the large kernel of the linear system of equations. The discretisation is based on the variationally equivalent method, as explained in subsection III.2.2.3. This allows the description based on the transfer operators N h and G h , so that the variationally equivalent operators can be defined using a projection approach, as detailed in equations (III.27) and (III.28). The discrete approximate inverse operators are denoted as

B Π ≈ X -1 h2 (III.56)
for operator (III.27) and

B G ≈ X -1 h3 (III.57)
for operator (III.28). The inverse of the original bilinear form X h1 is approximated by a smoothing operation

R h ≈ X -1
h1 .

(III.58)

The resulting algebraic multigrid preconditioner is based on a one-cycle reduction technique, which enables a fast application. The memory demand and setup time are reduced, since only one coarsening level is featured in the multigrid step.

III.4.1 Impact of coarsening type on convergence behaviour

The choice of coarsening operator has a large impact on the quality of the multigrid operator. The coarsening technique impacts the quality of the coarsening operator and, therefore, the convergence behaviour of the solution phase. Figures III.3 and III.4 show some typical combinations for the variationally equivalent operators. The discretisation is largely affected by the granularity of the partitioning.

A large number of partitions, which is equivalent to a high processor count, results in a high ratio of internal boundaries towards unique partitioning volumes. Therefore, the inter-process communication is very high for highly parallel test cases, as presented in figure III.4. In contrast, lower processor counts result in large partitioning volumes with relatively small processor boundary surfaces, as shown in figure III.3. Both figures contrast different combinations of coarsening methods for the variationally equivalent operators (III.27) and (III.28). The first abbreviation in the legend specifies the choice for operator (III.56), whereas the second label indicates the coarsening method for operator (III.57), which is related to the large null space of the linear system of equations.

For both figures it is visible that the convergence behaviour is largely impacted by the choice for operator B G that is used to reduce the error associated to the large null space. Even though operator B Π has a measurable impact on the convergence behaviour, it is shown that the correct description of the kernel is important for achieving a high error reduction rate.

The efficiency of the coarsening operators is mostly impacted by the surface area of the inter-processor boundary surfaces. For the lowly parallelised test case, shown in figure III.3, it can be noticed that the RGST, CLJP, FALG, as well as the, ECGC method show a comparable convergence behaviour. The RGST method naturally shows a very good convergence behaviour for this test case, since the setup structure is nearly sequential. The RGST method results in operator descriptions, with very high quality interior points. The discretisation quality only diminishes in the boundary regions. The lowly parallelised test case can, therefore, naturally be described by this method. The CLJP method has been developed to address the bad quality of the coarsening in the surface regions. Its use results in higher coarsening qualities for large numbers of partitionings with respect to the RGST method. Unfortunately, the interior coarsening quality is slightly reduced, when compared to the classical RGST technique. For the lowly parallelised test case on four processors, both methods show a near identical behaviour. The behaviour of the FALG scheme is naturally equivalent, because it is a combination of the RGST and CLJP methods. The ECGC scheme is used for very high processor counts and addresses the problem of very high parallel granularity. For the low partition count, depicted in figure III.3, it is visible that the behaviour is very close to the classical coarsening schemes. Its coarsening quality is nearly identical to the RGST, CLJP, as well as the FALG scheme. This is in stark contrast to the HMIS method, which results in weak convergence for all operator combinations. The use of the HMIS method results in coarsened operators with low operator complexities. It is used to reduce the memory demand and the cost of application, due to the reduction of coarse reconstruction points. The side effect is a large reduction in the convergence rate.

For the highly parallel test case, shown in figure III.4, the differences for each operator combination are stronger. The solver combination using the CLJP method shows the strongest error reduction, whereas the HMIS method leads to the weakest convergence. For this highly parallel test case it is noticeable that the RGST method results in an increase of iteration counts, when compared to the CLJP, FALG, as well as the ECGC techniques. This can be explained by the bad coarsening quality on the processor boundaries.

The results indicate that the most important element in the discretisation phase is the adequate representation of the variationally equivalent operator (III.28) that addresses the large null space of the linear system of equations. The impact of operator (III.27), which is based on the transfer operators for nodal to edge values is measurable, but clearly offset by the coarsening method for the former operator. Preconditioned relative residual vs. number of iterations The choice of coarsening method impacts both the setup and the solution phase. The setup phase is an initial cost that impacts the overall solution time and is an important aspect for problems with varying material properties, due to the need for reevaluating the entries in the linear system of equations. The solution phase is related to the application cost of the operator and the operator application count.

A high error reduction rate, which leads to a low iteration count, is beneficial to a reduction of the solution time.

Figures III.5 and III.6 show the setup and solution times for a 4 processor and a 16 processor test case. In both instances, it can be observed that the setup time is very much affected by the choice for the variationally equivalent operator B Π , which is much larger than operator B G . The operator that is related to the kernel of the linear system of equations has no practical influence on the setup time, due to its much smaller discretisation size.

For the 4 processor test case it can be seen in figures (III.5a-III.5e) that the HMIS method leads to faster setup times, compared to the CLJP, RGST, FALG and ECGC method. The RGST method is the slowest, due to the problematic treatment of inter-processor boundary nodes. The CLJP, FALG and ECGC method are nearly equivalent.

The solution time is visibly related to the coarsening operator B G , because it largely affects the convergence behaviour of the multigrid solver. The choice of the HMIS method for this operator increases the iteration count and, therefore, the solution time, as can be seen in figures (III.5f-III.5j). The strongest increase in solution time is visible in figure III.5g for the combination of RGST solver for operator B Π and HMIS method for operator B G . In contrast, it is very advantageous to employ the cost effective HMIS method for the discretisation of the larger operator, in combination with a strongly converging discretisation method for operator B G , as can be seen in figure III.5i.

The impact of domain partitioning can be seen in figures (III.6a-III.6e). The very ineffective RGST method, shown in figure III.6b, has a setup cost that is four times larger than the HMIS method shown in figure III.6d. The CLJP and FALG method, shown in figures III.6a and III.6c show a comparable behaviour. The advantageous use of both CLJP and RGST method in case of the FALG method efficiently reduces the setup time and minimises the discretisation problems of the RGST method in this highly partitioned test case. Figure III.6e demonstrates that the ECGC method is almost as efficient as the HMIS method.

It is again recognisable that the application cost has a large impact on the solution time. Figures show that the HMIS method should be avoided as a discretisation choice for operator B G , whereas it can be advantageous to use for operator B Π , as can be seen in figure III.6i. Figure III.6g demonstrates that it is highly disadvantageous to combine the RGST method and the HMIS method, due to the large iteration count and high application cost. The solution time for this combination is roughly twice as costly as for the more efficient combinations. It can be noticed that the ECGC, CLJP, FALG and the HMIS method are an effective choice for operator B Π . An advantageous choice of discretisation combinations, i.e. a strongly converging operator for the operator related to the null space of the linear system of equations and a low cost operator for the larger operator, can reduce the computational time by roughly 50% to around 40 s, when compared to the worst case combination.

R G S T , C L J P R G S T , R G S T R G S T , F A L G R G S T , H M I S R G S T , E C G C (g) 0 20 40 60 80 100 F A L G , C L J P F A L G , R G S T F A L G , F A L G F A L G , H M I S F A L G , E C G C (c) 0 50 100 150 200 250 F A L G , C L J P F A L G , R G S T F A L G , F A L G F A L G , H M I S F A L G , E C G C (h)
E C G C , C L J P E C G C , R G S T E C G C , F A L G E C G C , H M I S E C G C , E C G C (e) 0 50 100 150 200 250 E C G C , C L J P E C G C , R G S T E C G C , F A L G E C G C , H M I S E C G C , E C G C (j)
R G S T , C L J P R G S T , R G S T R G S T , F A L G R G S T , H M I S R G S T , E C G C (b) R G S T , C L J P R G S T , R G S T R G S T , F A L G R G S T , H M I S R G S T , E C G C (g) 0 20 40 60 80 F A L G , C L J P F A L G , R G S T F A L G , F A L G F A L G , H M I S F A L G , E C G C (c) F A L G , C L J P F A L G , R G S T F A L G , F A L G F A L G , H M I S F A L G , E C G C (h)
E C G C , C L J P E C G C , R G S T E C G C , F A L G E C G C , H M I S E C G C , E C G C (e) E C G C , C L J P E C G C , R G S T E C G C , F A L G E C G C , H M I S E C G C , E C G C (j)

III.4.3 Impact of coarsening threshold on convergence behaviour

The size of the coarsened structure is principally influenced by the coarsening strength threshold α t , as defined in equation (III.50). A high coarsening threshold parameter leads to a large drop of coarsening information, which results in a low operator complexity. The discrete operators generally have a smaller memory demand, but also lead to a reduction of the convergence rate. Figures III.7 and III.8 show the effects of using different combinations of operators with different operator complexities on the convergence behaviour, the setup, as well as the solution time for the solution stage on 4 and 16 processors. The coarsening algorithm that is used in this test case is a standard HMIS approach using a one-cycle AMS preconditioner. The first threshold in the legend denotes the strength threshold for operator B Π , whereas the second number indicates the strength threshold for operator B G . In both figures, the impact on the convergence behaviour, the setup time and the solver time is shown.

The results for the 4 processor test case, shown in figures (III.7a-III.7c) show again a strong dependency on the coarsening quality of the operator, which is related to the large null space of the linear system of equations. The coarsening quality of the other operator has only minor influence. For operator values of α t = 0.25 and α t = 0.5 the convergence behaviour is nearly identical, which is also reflected in the setup times and the solution times, as shown in figures III.7d and III.7e. It is visible that a very high strength threshold, i.e. a low operator complexity, for operator B G leads to a large reduction of the convergence rate, which is shown in the greater solver times. The advantage of a bigger strength threshold is a reduction in memory demand. Figure III.7d demonstrates that a reduction of the operator complexity for operator B Π measurably reduces the setup cost for α t = 0.9.

The same behaviour is demonstrated for the highly partitioned test case on 16 processors in figure show a strong dependency of the convergence rate and the good quality reproduction of the operator related to the large null space of the linear system of equations. The qualitative reduction of setup time in comparison to the lowly parallelised test case is shown in figure III.8d, which shows that the most efficient setup is possible, if operator B Π is discretised using a low operator complexity. The solution times are reduced, but only by approximately 35%. It is again visible that a large increase of the strength threshold for operator B G leads to a huge increase in solver time.

Both test cases demonstrate that the strength threshold can safely be increased to α t = 0.5 without large negative effects on the convergence rate, setup or solution time. The advantage of reducing the operator complexity is a reduction of memory demand for the storage of the discrete operator structures. It is noticeable that the good coarsening quality of the operator dealing with the error reduction associated to the large null space is again the most important aspect of a good error reduction rate, which is also demonstrated in subsections III.4.1 and III.4.2. 10 -1 10 0 10 1 10 2 10 3 0.9,0.25 0.9,0.5 0.9,0.9 0 .2 5 ,0 .2 5 0 .2 5 ,0 .5 0 .2 5 ,0 .9 0 .5 ,0 .2 5 0 .5 ,0 .5 0 .5 ,0 .9 0 .9 ,0 .2 5 0 .9 ,0 .5 0 .9 ,0 .9 (d) solver setup time in s 0 50 100 150 200 250 0 .2 5 ,0 .2 5 0 .2 5 ,0 .5 0 .2 5 ,0 .9 0 .5 ,0 .2 5 0 .5 ,0 .5 0 .5 ,0 .9 0 .9 ,0 .2 5 0 .9 ,0 .5 0 .9 ,0 .9 10 3 0.9,0.25 0.9,0.5 0.9,0.9 0 .2 5 ,0 .2 5 0 .2 5 ,0 .5 0 .2 5 ,0 .9 0 .5 ,0 .2 5 0 .5 ,0 .5 0 .5 ,0 .9 0 .9 ,0 .2 5 0 .9 ,0 .5 0 .9 ,0 .9 2 5 ,0 .5 0 .2 5 ,0 .9 0 .5 ,0 .2 5 0 .5 ,0 .5 0 .5 ,0 .9 0 .9 ,0 .2 5 0 .9 ,0 .5 0 .9 ,0 .9 The typical application of the AMS method involves a symmetric application of smoothing and multigrid operators. The normal configuration involves an error smoothing operation R h , the error reduction for the projected vector components B Π , defined in equation (III.56) and the reduction of the residual error for the large null space of the linear equation system, denoted by B G , as defined in equation (III.57).

The major focus of the standard AMS method is the reduction of the error due to the large null space of the linear system of equations. Its focus lies, therefore, on the application of operator B G . This is normally achieved in a multiplicative manner, but can also be done using an additive form, here denoted as AMS-A. An alternative implementation of this complete preconditioning operator is denoted as AMS-Π, which focusses on operator B Π . In addition, it might be advantageous to smooth out the error after each application of the coarsened operator using additional smoothing cycles R h . Table III.1 gives an overview of different operator orderings that can be easily implemented based on the coarsened operators. III.1: Different operator orderings for the AMS method Figure III.9 shows the impact of the different operator orderings on the convergence behaviour for a preconditioning operator, which is based on an HMIS coarsening method. It is visible that the additive version of the standard preconditioner results in a large reduction of the convergence rate. The error is very slowly reduced, which leads to a large increase in computational effort. The additive form leads to an increase in computational time by a factor of approximately 2. The AMS-Π preconditioner converges slower than the standard AMS solver. In the above-mentioned subsections it is demonstrated that a large error reduction is only possible, if the preconditioning operator focusses on the large kernel of the linear system of equations. It can be seen that the overweighting of operator B Π with respect to B G leads to an increase in computational time by roughly 50%. In addition, it can be noted that additional smoothing operations do not lead to an increase in the error reduction rate. In fact, the computational time is increased, due to an increase in computing cost, without a visible increase of the convergence rate, as is shown for the AMS-GS ordering type.

type operator ordering application type

AMS R h , B G , B Π , B G , R h multiplicative AMS-A R h , B G , B Π , B G , R h additive AMS-Π R h , B Π , B G , B Π , R h multiplicative AMS-GS R h , B G , R h , B Π , R h , B G , R h multiplicative Table

III.4.5 Impact of material parameter distribution on convergence behaviour

A typical steel work piece can show a highly temperature dependant behaviour.

The relative magnetic permeability is greatest for the initial temperature and decreases for rising temperature values. The changes of the material parameters are reflected in the values of the linear system of equations, which are weighted differently according to the material parameters at the integration points inside each finite element. In addition, large material jumps are possible, due to the constant permeability of the surrounding air, which is fixed at µ r = 1, whereas the relative magnetic permeability of steel can be several magnitudes higher. The large material discontinuities can lead to instabilities in the solution phase, due to an impact on the convergence quality. In order to model the impact of temperature dependency, a reasonable boundary test case is shown. Initially, it is assumed that the crankshaft possesses a constant homogeneous relative magnetic permeability of µ r = 10. The computation is advanced by one predetermined time step and the resulting temperature distribution is utilised to evolve the relative magnetic permeability. Here, the function µ r = f (B) is defined using a predetermined range beginning by µ r = 1 and increasing to the initial relative magnetic permeability µ r = 10. It is assumed that the maximum temperature is equal to the specific Curie temperature of the material (cf. subsection I.5.2). The temperature range is further subdivided into ten equal parts, which are each associated to a corresponding relative magnetic per-meability. This stepwise association describes a reasonable boundary test case for temperature dependant material evolution models. In practice, the temperature distribution, as well as the resulting relative magnetic permeability, are smooth, due to the smooth transitioning of the temperature distribution. Here, the described test case leads to strong discontinuities inside the work piece and large material jumps. These jumps may lead to instabilities that can be assumed to model a worst case scenario.

Figure III.10 shows the discontinuous material distribution for a cut of the the supporting sides of the crank shaft domain. It is visible that the material evolution leads to large material jumps, including large differences near the problematic boundary edges. The highest material jumps are, therefore, near the parts of the computational domain that are most difficult to treat numerically.

Figure III.11 shows another view for the heat treated region. It is again visible that the chosen material evolution model leads to large discontinuities with a complicated structure. The step wise layering leads to large material jumps inside of the crankshaft domain and, therefore, directly affects many neighbouring entries in the linear system of equations. .12 gives results for the convergence behaviour using a standard HMIS based AMS preconditioned solver. It is visible that an increase in the relative magnetic permeability and, thus, a difference in weighting of the entries in the linear system of equations, leads to a reduction of the convergence rate. Yet, a very large increase of the material jumps by a magnitude of 3 for µ r = 10 3 leads only to a slight increase in computational effort by a factor of approximately less than 2. The test case that is based on the highly discontinuous material In effect, the many layers of discontinuities do not lead to a reduction of the convergence rate compared to the test case with constant relative magnetic permeability of µ r = 10. It can be observed that the AMS preconditioned solver is not destabilised and converges monotonically. For every test case, the solution method is stable and consistent.

IV Modelling inductor motion

Français:

Le chapitre suivant décrit une nouvelle approche pour le mouvement de l'inducteur dans une configuration Lagrangienne en utilisant une transformation rigide de la surface de l'inducteur.

English:

The following chapter describes a novel approach for the inductor movement in a Lagrangian setting using a local remeshing technique that uses a rigid transformation of the initial inductor surface. The identification function is based on a discrete level set function, which is created using the surface primitives of the discrete inductor domain. The resulting algorithm is local, efficient and can be parallelised.

In the beginning, the problem is outlined. Afterwards, an overview of existing methods in the context of induction heating and other applications is given. Classical approaches using the Arbitrary-Lagrangian-Eulerian (ALE) method, as well as discrete algorithms are presented. The advantage of a local level set approach is shown. Following is the description of the proposed algorithm. A two-dimensional example domain has been chosen to illustrate the level set identification technique. Afterwards, results are presented for a test case: the rotation of an inductor around a crankshaft. It is shown that the initial mesh quality is preserved in subsequent time steps.

IV.1 Problem statement

It is mentioned in section I.3 that the industrial application of the heat treatment process of the automotive crankshaft involves a relative motion of inductor and crankshaft during the procedure. The crankshaft is placed in a mounting device, visible in figure IV.1 that allows a relative rotation during the heating stage. Afterwards it can be used to immerse the workpiece in a cooling liquid, like water or oil or it can act as a support for cooling showers. The relative movement is included into the heat treatment procedure for several reasons. The eddy currents are created by the rapidly changing magnetic field, circulating in the surface of the crankshaft domain. The geometric centre of this loop will not feature any measurable current, due to current cancellation effects. A heating in these regions can only be accomplished, either by conduction or by a relative movement of the inductor, so that the region is shifted to a location without current cancellation effects.

The skin effect results in a distribution of eddy currents near the surface of the domain. An approximation for the penetration depth of a homogeneous domain has been presented in section I.5. The exponential decaying model is presented in equation (I.6). It can be seen that the penetration depth depends inversely on the applied frequency, i.e. a high frequency leads to a small penetration depth. The higher the frequency, the higher any heating pattern error will be, if the inductor is slightly misaligned, with respect to the conductor. In addition, it must be noted that neither the inductor geometry nor the workpiece is fully symmetric. A relative movement during the heat treatment can help to smooth out the applied heating energy. An industrial application using relative movement can increase the reliability of the induction heating process.

The overall goal of an induction heat treatment of the crankshaft is the generation of a layer of austenite near the contact surface, where the bearing will be placed. This austenitic phase will be transformed into the martensitic state in a subsequent quenching procedure. The relative motion ensures that the heat can be applied evenly, even though the geometry is relatively complex. E.g. the automot-ive crankshaft includes problematic features like holes and sharp corners, which can lead to either under-or over-heating during the induction process. The heat treatment process is, in essence, an interaction of the electromagnetic induction heating, which is the major driving force for the increase of the temperature, as well as heat conduction and relative motion, which both allow the heat to distribute evenly and towards regions that are not directly affected by strongly induced eddy currents. For these cases, heat conduction is advantageous with respect to an application of heat in these regions, even though it negatively affects the centre of the workpiece, i.e. in general, conduction is avoided, since it leads to an increase in temperature in the interior of the workpiece. Therefore, the induction heating process is achieved as quickly as possible.

The introduction of relative movements in numerical simulations is difficult, since it complicates the numerical discretisation. Existing numerical models like the remeshing applications for forming processes or moving domain problems in fluid dynamics are highly problem dependent and can not be used directly for the induction heating process in an efficient manner. This is due to assumptions made in these models that might not be applicable to the induction process or features that might lead to a severe increase in computational complexity, such that the practical application is not possible in this context. Section IV.2 gives an overview over different concepts of general models for relative movements in numerical simulations, as well as in the context of induction heating. Afterwards, section IV.3 proposes a new model for the description of inductor motion in the context of electromagnetic induction heating that is efficient and parallelisable. It is followed by section IV.4, which presents results for the application of inductor motion in the context of a global finite element analysis including the geometry of a sample industrial crankshaft.

IV.2 Existing methods

Existing methods for the description of movements of computational domains can be divided into several sub-categories. First, it is illustrative to sort them by the approach that is used to describe the movement. Generally, there are two configurations that are applied in classical mechanics. The first description is the Lagrangian approach, which attaches quantities of interest to a particle, which moves in time, i.e. any quantity like the temperature or the density is attached to a physical object that might change its shape and position through time. In contrast, the Eulerian approach fixes a certain control volume at a spatial location. The quantities of interest are allowed to pass through at any given time step. The first approach associates the quantities to objects, whereas the second approach associates quantities to spatial locations.

The Lagrangian approach is advantageous, since it allows to consistently track the outlines of the real geometries. Unfortunately, it is very difficult to include movements in such a description, since these generate large distortions in the grid. The Eulerian approach observes the physical quantities, moving in a fixed grid, therefore, the mesh does not get distorted. The disadvantage of this approach is that the boundary surface of the approximated object might not align with the boundary of the discretisation, which leads to interface approximation errors.

A generalisation of the above-mentioned methods is the Arbitrary Lagrangian Eulerian (ALE) approach that has been introduced to utilise the advantages of both approaches, while minimising the shortcomings. A general introduction to this method can be found in [START_REF] Donea | Arbitrary Lagrangian-Eulerian Methods[END_REF]. In the ALE approach neither the material (Lagrangian) nor the spatial (Eulerian) configuration is used solely for the description of the evolution of the physical quantities. Instead, a third intermediate reference configuration is employed. This reference configuration is used to track the computational grid during the evolution of the boundaries. In this approach, small movements are represented using an update of the grid itself, stemming from a Lagrangian point of view, whereas large deformations, which would lead to a large distortion of the discretisation, are represented using a Eulerian approach. The ALE method has been used extensively in moving boundary applications, e.g. [START_REF] Nguyen | An arbitrary Lagrangian-Eulerian discontinuous Galerkin method for simulations of flows over variable geometries[END_REF] and [START_REF] Persson | The Numerical Simulation of Flapping Wings at Low Reynolds Numbers[END_REF], which use ALE to describe flow around deformable domains or [START_REF] Feldman | Fluids in Deforming Meshes[END_REF] that uses an ALE approach in a mostly Eulerian fluid flow formulation, allowing mesh deformations using the semi-Lagrangian surface tracking technique of [START_REF] Bargteil | A Semi-Lagrangian Contouring Method for Fluid Simulation[END_REF] to incorporate advection velocities in a deformed mesh.

Unfortunately, the ALE method can pose problems with respect to approximation errors, due to interpolation errors if large deformations are computed. E.g. [START_REF] Wall | Large Deformation Fluid-Structure Interaction -Advances in ALE Methods and New Fixed Grid Approaches[END_REF]] mentions for fluid structure interaction problems that remeshing must be done, if the underlying discretisation is severely distorted, which can arise for certain problems, stating "Such scenarios include problems with topology changes, e.g. when a fluid penetrates an opening crack or when a valve closes, or simply if the structure moves too far in the domain or rotates."

In addition to the categorisation in material or spatial tracking approaches, it can be useful to divide the employed methods in further sub-categories based on the choice of discretisation. The approaches mentioned above, i.e. Lagrangian, Eulerian, and the ALE technique are valid in a continuous, as well as a discrete context. These approaches can be used for various discretisation methods. Examples are methods utilising a weak formulation in a Galerkin finite element context, finite difference approximations or boundary element approaches. E.g. the boundary element approach can be used in a Lagrangian context, mostly in a mixed formulation due to the necessity of fully discretising the conductive domains. Representations for this approach with respect to an application of moving domains in electromagnetics can be found in [START_REF] Frangi | Coupled Fast Multipole Method -Finite Element Method for the analysis of magneto-mechanical problems[END_REF][START_REF] Frangi | Magneto-mechanical simulations by a coupled fast multipole method-finite element method and multigrid solvers[END_REF], whereas general mixed boundary element formulations can be found in [START_REF] Alotto | Three-dimensional eddy current analysis in unbounded domains by a DEM-BEM formulation[END_REF][START_REF] Camaño | Analysis of a FEM -BEM model posed on the conducting domain for the time-dependent eddy current problem[END_REF][START_REF] Meddahi | An H-based FEM-BEM formulation for a time dependent eddy current problem[END_REF][START_REF] Rodríguez | A FEM-BEM approach for electromagnetostatics and time-harmonic eddy-current problems[END_REF]. [START_REF] Frangi | Coupled Fast Multipole Method -Finite Element Method for the analysis of magneto-mechanical problems[END_REF][START_REF] Frangi | Magneto-mechanical simulations by a coupled fast multipole method-finite element method and multigrid solvers[END_REF] describe the mechanical problem of an industrial relay using a mixed FEM/BEM formulation, allowing relative movement without the need of remeshing. [START_REF] Alotto | Three-dimensional eddy current analysis in unbounded domains by a DEM-BEM formulation[END_REF] solves a general eddy current system without imposing movements using a discrete cell element model for the conducting regions, coupled using boundary elements. A movement could easily be incorporated, which is also possible for the general mixed-models proposed in [START_REF] Camaño | Analysis of a FEM -BEM model posed on the conducting domain for the time-dependent eddy current problem[END_REF][START_REF] Meddahi | An H-based FEM-BEM formulation for a time dependent eddy current problem[END_REF][START_REF] Rodríguez | A FEM-BEM approach for electromagnetostatics and time-harmonic eddy-current problems[END_REF]. The advantage of these approaches is that a relative movement, as well as the vanishing magnetic field condition at a certain distance from the conducting objects, can easily be taken into account. The disadvantage is the need to solve a coupled linear system of equations that is difficult to handle numerically (cf. chapter III).

Another discretisation based approach has been presented in [START_REF] Alotto | Efficient Use of the Local Discontinuous Galerkin Method for Meshes Sliding on a Circular Boundary[END_REF], which introduces a local discontinuous Galerkin method for the computation of rotating electric machines in a global embedding domain. The method is based on introducing a slip-surface that allows hanging nodes, such that a relative circular movement of a rotor with respect to a stator can be included in the discretisation. The discontinuous Galerkin method utilises a finite element description in such a way that the solutions are allowed to have discontinuous values across element boundaries. The inter-element connection is established by incorporating a numerical flux, arising in the weak-form (cf. section II.2.2.3) that approximates the real inter element flux. The freedom of choice for this numerical flux makes the discontinuous Galerkin method very applicable to hyperbolic conservation problems. A general introduction to this method can be found in [START_REF] Hesthaven | Nodal discontinuous Galerkin methods: algorithms, analysis, and applications[END_REF] or [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF]]. An example for the discontinuous Galerkin method, applied to problems including second order derivatives in space, has been presented in [START_REF] Cockburn | The Local Discontinuous Galerkin Method for Time-Dependent Convection-Diffusion Systems[END_REF]] (cf. [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF] for a general discussion of the discontinuous finite element method for elliptic problems). This flux-based coupling enables the introduction of hanging nodes in the finite element mesh, as long as the inter-element fluxes can be prescribed in a consistent manner. [START_REF] Alotto | Efficient Use of the Local Discontinuous Galerkin Method for Meshes Sliding on a Circular Boundary[END_REF] suggests this method can be used to deal with the moving mesh problem in electromagnetics. Unfortunately, the hanging node approach is only feasible for domains including a certain symmetry, making the movement of complex domains problematic, e.g. the example in [START_REF] Alotto | Efficient Use of the Local Discontinuous Galerkin Method for Meshes Sliding on a Circular Boundary[END_REF] features a circular shaped rotor. In addition, the introduction of discontinuous elements increases the degrees of freedom of the numerical system. Furthermore, the hanging nodes can result in solutions featuring spurious oscillations on irregular meshes, requiring additional correction methods, as shown by [START_REF] Buffa | The Mortar-Discontinuous Galerkin Method for the 2D Maxwell Eigenproblem[END_REF].

Another discretisation technique for incorporating moving domains in the eddy current formulation is the dual mesh technique proposed by [START_REF] Flemisch | Coupling scalar and vector potentials on nonmatching grids for eddy currents in a moving conductor[END_REF]. It incorporates the movement of the conductive domains by using non matching overlapping triangulations using a mortar element technique. Unfortunately, the projection operator used in this context is based on a Lagrange multiplier technique, which is numerically difficult to handle.

An example for a Eulerian approach can be found in [START_REF] Bay | A numerical model for induction heating processes coupling electromagnetism and thermomechanics[END_REF]] that introduces the movement of the inductor, by utilising an inter-element transportation method of the inductor material in a fixed mesh. This standard Eulerian technique is problematic in the the context of induction heating, since it is very important to track the correct boundaries of each object. Otherwise, the discretisation might not conform to the physical problem, i.e. the surface effect could lead to a large deviation of the induced eddy currents inside the domain, leading to further approximation errors, which might amplify in subsequent time steps.

It follows that the simplest description in the context of induction heating is a standard Lagrangian approach. This allows a tracking of the initial interface, ensuring that the physical effects can be correctly reproduced. Unfortunately, in practice it is problematic to include the moving boundary problem in such a global Lagrangian description. In essence, the finite element mesh has to be adapted at every time step of the numerical computation, in order to reflect the correct relative placement of each respective physical object. A standard approach in this context is the mesh adaptation of an existing finite element description, based on a remeshing indicator.

An example for a remeshing technique of this type can be found in [START_REF] Boussetta | Adaptive remeshing based on a posteriori error estimation for forging simulation[END_REF]] that defines a remeshing technique based on a Zienkiewicz-Zhu errorestimator [START_REF] Zienkiewicz | A Simple Error Estimator and Adapative Procedure for Practical Engineering Analysis[END_REF]. It defines an adaptive remeshing technique for an existing mesh that undergoes a large deformation, imposing a local mesh improvement procedure during the process. The examples provided in [START_REF] Boussetta | Adaptive remeshing based on a posteriori error estimation for forging simulation[END_REF] include workpieces that are deformed by external tools. The boundary conditions are provided using surface-to-surface contact methods. The material domain is described using a full Lagrangian approach, so that the material itself is tracked using the elements of the finite element description. The model proposed in [START_REF] Boussetta | Adaptive remeshing based on a posteriori error estimation for forging simulation[END_REF]] deals, therefore, with separate domains that are not embedded in a global finite element domain, but can be treated independently, due to the contact condition. In contrast, the induction heating problem in this work deals with a configuration of workpiece and surrounding inductor that is embedded in a surrounding air domain. A remeshing method, like the one proposed in [START_REF] Boussetta | Adaptive remeshing based on a posteriori error estimation for forging simulation[END_REF]] could be used to describe this problem, even though a fully global remeshing operation might be numerically difficult to apply. In addition, mesh coarsening and refining operations can lead to interpolation and approximation errors that lead to model errors. A remeshing process that is based on error estimators can not guarantee the consistency of the discretised surface during the computation, i.e. for the problem of induction heat treatment the remeshing might lead to a discretisation error of the inductor or the crankshaft surface, which can have a severe impact on the quality and form of the computed magnetic field.

It is demonstrated in chapter III that the electromagnetic computation of the induction heating problem can be efficiently solved, if a global surrounding air domain approach is used. This method leads to large, but sparse symmetric matrices that can be solved using very efficient solution techniques. Therefore, a very advantageous approach should strive to include the relative movement of inductors in such a manner that every subsequent time step is essentially comparable to the first time step, in terms of mesh discretisation and mesh quality. In addition, a prospective method should not necessitate the introduction of additional unknowns in the final electromagnetic problem, which might destroy the good convergence properties or the applicability of the efficient solvers that are presented in chapter III.

In the following section IV.3 a mesh identification and mesh readaptation algorithm is proposed that can be used to describe the inductor motion in a global finite element mesh using discrete level set functions. The fundamental observation is that the interfaces of the objects do not change their shape or come into contact during the heat treatment. Therefore, any description of motion can be described using rigid body transformations in a Lagrangian setting in such a way that the initial mesh description reflects the correct boundaries throughout the full computation. The algorithm is based on the discrete representation, i.e nodes, edges and faces of the discrete elements. The practical application is based on geometrical searching algorithms that can be implemented in an efficient manner and are local in nature, so that they can be treated in parallel. The resulting finite element mesh has the same quality as the initial mesh in every subsequent time step of the numerical computation.

IV.3 Inductor motion using discrete level set functions

The proposed algorithm is presented using a two-dimensional example domain, depicted in figure IV.2a. It consists of an embedding model domain Ω that fully encloses an inductor Ω tn I . The method is indifferent to the number of enclosed domains, e.g. inductor, crankshaft and flux concentrator. Therefore, it suffices to demonstrate the algorithm using a simple configuration with only one conducting element. A sample triangulation of the initial computational domain Ω using triangular finite elements is visible in figure IV.2b. The full computational domain is given as T , whereas the inductor domain is specified as T tn I . The index t n specifies that this domain is an initial time step. The initial geometry serves as basis for the transformed domain T t n+1 I that will be used in the next time step of the computation. The following time step, depicted as t n+1 is the next discrete time increment for the imposed movement. The key observation for induction heat treatment is that the surfaces of the involved objects can be assumed to be rigid, so that the initial domains can change their position and orientation but not their shape. It follows that the transformed inductor domain Ω t n+1 I can be described using a rigid body transformation, as

Ω t n+1 I = Ψ(Ω tn I , t n+1 ), (IV.1)
which can directly be applied to its discrete form T t n+1 I

, such that

T t n+1 I = Ψ(T tn I , t n+1 ). (IV.2)
Another observation is that neither contact nor inter-object penetrations can occur during the heat treatment. In addition, each rigid object moves in the embedding air, which is not affected by the electric computation, i.e. neither a temperature change nor the applied magnetic loading are assumed to change its relative magnetic permeability or its conductive properties. This invariant property of the air domain leads to the conclusion that the movement of the inductor passes entirely through a domain of which the properties are known in advance for each subsequent computational step. The main idea of the proposed algorithm is that it is theoretically possible to fully remove the embedding air domain, followed by a connecting remeshing operation. This idea originates from observations regarding the creation of the initial finite element discretisation. In the beginning, several conducting domains are placed in a surrounding container object. Afterwards the air-gap is filled with finite elements, in order to create a conforming finite element mesh for the electromagnetic computation. The conclusion is that it would be possible to initially transform the inductor in the first time step, before the computational mesh has been processed, followed by the connecting meshing step. The resulting finite element mesh would qualitatively be comparable to the nontransformed initial time step.

The disadvantage of this approach is that a full meshing of the embedding domain is computationally intensive. Furthermore, the meshing time, as well as the quality of the resulting mesh are undetermined for arbitrary transformations. In practice, the initial finite element mesh is prescribed by meshing algorithms, like the ones presented in [START_REF] Peraire | Unstructured mesh generation including directional refinement for aerodynamic flow simulation[END_REF] or in [START_REF] Gruau | 3D tetrahedral , unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric[END_REF], which are assisted by metrics, specifying a prescribed element size. These bounding geometries have to be described manually, in order to ensure the correct initial meshing of complicated features. It is furthermore useful to initially prescribe the element size in areas that are important for the numerical analysis. E.g. a crankshaft that is heat treated by one rotating inductor can be meshed, such that areas with strong gradients are finely meshed, whereas not strongly affected features in a certain distance to the inductor are meshed in a coarse manner, in order to save memory and computational time.

The proposed algorithm rests on the assumption that an initial finite element discretisation is valid, such that any element in a certain distance to the initial geometries will have been discretised using element sizes that are consistent with the physical model. Let meas define a measure of a geometry, such that the intersection of the transformed geometry Ω t n+1 I and initial geometry Ω tn I can be used to bound the relative change of the permissible range of the rigid body transformation. The prescribed bound is defined by meas(Ω

t n+1 I ∩ Ω tn I ) meas(Ω tn I ) ≥ Z, (IV.3)
where Z is a user supplied constant

0 ≤ Z ≤ 1, (IV.4)
which bounds the relative permissible change of the transformable geometry. If Z is sufficiently large, the transformation will move through a layer of the embedding air-domain that will possess a comparable element size. The goal is to identify this layer of elements, which can be removed, due to the assumption of the above-mentioned invariant properties of the constitutive elements. This surface element layer defines a boundary that can be used akin: to the initial boundary domain to recreate a conforming finite element mesh. The key difference is that the finite element description can be assumed to include enough information to facilitate a remeshing step, without the need of manual interaction, like in the above-mentioned global remeshing case. In practice, the parameter Z will depend on the rigid body transformation, as well as the initial finite element mesh. A choice of Z = 0.8 would indicate that the permissible change is valid, if the transformation of the boundary will result in a transformed geometry that coincides with 80 % of the volume of the initial step. The condition, described in equation (IV.3) can be used throughout the computation, for validating the imposed time step. If the condition indicates that the transformation exceeds the permissible range, the time step can simply be reduced. This does not change the outcome of the computation, since the identification algorithm is based on the discrete form of the correct boundaries. A smaller time step that passes condition (IV.3) will necessarily result in a consistent finite element description.

The removable subset of discrete elements will be denoted as T t n+1

E R . The index denotes that these elements are related to the embedding air domain. The transformed inductor domain T

t n+1 I
is embedded in this air domain. The removable subset can be described using an identification function L (x) that describes, whether a discrete geometrical feature, represented by the geometrical point x, is affected by the transformation described using equation (IV.2). The removable subset can, therefore, be defined as

T t n+1 E R = (T tn I ∪ T t n+1 L ) \ T t n+1 I , (IV.5)
using the discrete subset of finite elements T t n+1 L

, which have been defined using L (x). It can be observed that the set of affected finite elements is only a small subset of all possible finite elements in the computational domain T . In practice, a further reduction technique can be applied to the embedding elements, such that

T t n+1 E R ⊆ T t n+1 E * ⊆ T tn E . (IV.6)
Figure IV.3a depicts the initial global finite element mesh, with an overlay of the transformed inductor domain T t n+1 I in red. The elements near the border, as well as the bottom range are unaffected by the transformation. In practice, a reduction of the searchable finite element space can be achieved by applying an embedding bounding geometry Ω . In practice, this operation can be achieved by an update on the vertices of the finite element mesh. The local element incidence information is not influenced by this transformation algorithm. Afterwards, the remaining mesh can be closed to create a conforming finite element mesh. The recreated mesh T

t n+1 R is depicted in violet in figure IV.5. It is a subset of the new embedding domain T t n+1 E
that can serve as a base for subsequent time steps of the inductor movement. The identification algorithm has led to a local remeshing strategy that affects only a small subset of finite elements of the computational domain. The characteristic element size of the affected elements should be in relation to the velocity and time step, which is also assured by the conformity condition described in equation (IV.3).

The mesh identification function L (x), which identifies the affected elements T t n+1 L , can be described independently of the local adaptation algorithm. In the context of induction heating with rotating inductors, it is advantageous to base this function on the initial transformed geometry, due to the assumption of object rigidity. In the following, this function is identified as a discrete level set function that is based on the geometrical features of the discrete inductor surface. In a twodimensional setting it consists of vertices and edges, whereas it consists of vertices, edges and faces in three dimensions.

The function can be described using the relative distances of these geometrical base features, with subsequent ordering according to the minimal distance. The resulting algorithm is simple and can be applied in parallel. Since it is based on the discrete version of the geometries, any identification data will be accurate up to the floating point accuracy. The mesh quality is not affected by the number of time steps, i.e. since the accuracy of the finite element mesh in subsequent time steps is always comparable to the initial time step, due to the restrictions of the permissible time step size and velocity, the accuracy remains good for all subsequent time steps.

A level set function is a spatial function that implicitly defines a surface. The function maps a geometrical point x in the scalar space of real functions R, defining whether it lies inside, outside or on the specified surface boundary. For simple geometries it can be stated explicitly. E.g. a circle in the two dimensional space can be specified by its center x c = [x c1 , x c2 ], as well as its radius r, such that L (x) can be defined for any point

x = [x 1 , x 2 ] as L (x) =  (x c1 -x 1 ) 2 + (x c2 -x 2 ) 2 -r. (IV.7)
Further implicit level set distance functions for simple shapes like rectangles, polygons, as well as boolean operations on these can be found in [START_REF] Persson | A simple mesh generator in matlab[END_REF] which proposes a mesh generation technique based on implicit level set distance functions. These types of level set distance functions can directly be used to identify the set T t n+1 E R . Unfortunately, it is difficult to describe a three-dimensional object like an inductor using polygonal shape functions. In general, level set functions are used in a Eulerian setting to describe the surface evolution. Such techniques are commonly employed in fluid dynamics, e.g. [START_REF] Barton | A conservative level-set based method for compressible solid/fluid problems on fixed grids[END_REF] or [START_REF] Marchandise | A stabilized finite element method using a discontinuous level set approach for the computation of bubble dynamics[END_REF]. The usual practice consists in computing an initial level set distance function L (x) in the setup phase. Afterwards, the level set function is transported, using a velocity field v. The general level set advection formula (see, e.g. [START_REF] Osher | Fronts Propagating with Curvature-Dependent Speed : Algorithms Based on Hamilton-Jacobi Formulations[END_REF] or [START_REF] Osher | Level Set Methods and Dynamic implicit Surfaces[END_REF]) relates the transient evolution of the level set surface to its gradient, such that

∂L (x) ∂t + v∇L (x) = 0. (IV.8)
This hyperbolic conservation equation is solved in subsequent time steps to compute the advective transport of the initial surface. It is numerically difficult to handle, thus requiring more elaborate solution techniques, like the streamlineupwind-Petrov Galerkin finite element methods that deal with discontinuities and shocks by introducing artificial diffusion in the weak form (see, e.g. [START_REF] Brooks | Streamline Upwind/Petrov-Galerkin Formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-Sokes Eequations[END_REF]) or discontinuous Galerkin finite element methods (see, e.g. [START_REF] Cockburn | TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Conservation Laws II : General Framework[END_REF]).

In this work it is proposed to describe a discrete version of the level set function based on the surface primitives of the initial inductor domain. This is possible, due to the assumption of rigidity. The level set function will, therefore, in the following be described by its discrete equivalent L (x, T

t n+1 I
) that is defined by

L (x, T t n+1 I ) = d(x, ∂T t n+1 I )sgn(n xs • ∆x).
(IV.9)

The symbol T

t n+1 I
denotes the set of finite elements that identify the discrete version of the inductor, whereas ∂T t n+1 I specifies the surface. In a three dimensional setting, these surface elements consist of triangles, if tetrahedral elements are chosen for the finite element discretisation. The first part of equation (IV.9) defines the infimum distance of the geometrical point x to the discrete surface ∂T t n+1 I , which can be defined as d(x, ∂T

t n+1 I ) = inf ∥ x -x s ∥, ∀s ∈ ∂T t n+1 I .
(IV.10)

The second part of equation (IV.9) specifies whether the geometrical point x is inside, outside or on the discrete surface. It uses the surface normal n xs of the closest geometrical feature, defined by point x s , of the discrete surface ∂T

t n+1 I
for the point x. The signum function is used to measure the sign of the scalar product between the normal vector n xs and the difference vector ∆x = xx s .

The general application of this discrete level set function, in order to identify the element subset T t n+1 I , does not necessitate the explicit knowledge of the distance itself. The definition of equation (IV.9) includes the scalar distance, since it is a by-product of the algorithm for searching the closest geometrical feature. Its value can be used in a remeshing step to monitor or impose local mesh sizes.

The discrete inductor surface is not C 1 -continuous. Therefore, there is no analytically defined normal at the vertices and edges of the discrete surface. As a remedy, a pseudo normal has to be implemented that approximates the continuous normal. E.g. [START_REF] Gouraud | Continuous Shading of Curved Surfaces[END_REF] proposes an average weighting of the normal n G xs , stating that "This normal could be computed as, for example, the average of the normals to each polygon associated with this particular vertex", such that

n G xs =  n i=1 n f i ∥  n i=1 n f i ∥ , (IV.11)
for the n connecting faces to the vertex point x s , each defining a surface normal n f i . Unfortunately, this simple averaging is problematic, since it equally weights each of the connecting face normals. Consider the examples in figures IV.7a and IV.7b. Equation (IV.11) defines the same surface normal for the vertex x s for the example presented in figure IV.7a as for the example in figure IV.7b. This leads to inconsistencies for complex finite element discretisations. The third tetrahedron it is out of view.

[ [START_REF] Baerentzen | Signed Distance Computation Using the Angle Weighted Pseudonormal[END_REF] proves that the angle weighted pseudo normal n α xs , provided in ( [START_REF] Séquin | Procedural Spline Interpolation in UNICUBIX[END_REF], p. 15), can instead be used for the sign computation in equation (IV.9). Comparisons to alternative pseudo normals, like the method of [START_REF] Glassner | Computing Surface Normals for 3D Models[END_REF] that creates a vertex normal based on a least squares minimisation technique or [START_REF] Huang | A Complete Distance Field Representation 1[END_REF] that defines the surface normal based on the largest inner product between the distance vector ∆x and the incident face normals, can be found in section 2.1 of [START_REF] Baerentzen | Signed Distance Computation Using the Angle Weighted Pseudonormal[END_REF].

The angle-weighted pseudo normal is a weighted version of the above-mentioned Gouroud normal, defined in equation (IV.11). Let α i be the incidence angle of the face, in relation to the vertex point. E.g. for the above-mentioned example in figures IV.7a and IV.7b the incident angle for the left face is α 1 , whereas it is α 2 for the right face. It follows that the angle-weighted pseudo normal is defined as the weighted sum

n α xs =  n i=1 α i n f i ∥  n i=1 α i n f i ∥ . (IV.
12)

The influence of the vanishing angle α 1 in figure IV.7b can therefore be accounted for. [START_REF] Baerentzen | Signed Distance Computation Using the Angle Weighted Pseudonormal[END_REF] shows that the angle weighted pseudo normal can be defined for vertices and edges, observing that it is a generalisation of the face normal, so that it can be used to distinguish whether a point lies inside or outside of a closed surface. Specifically, the following inequalities hold sgn(n xs • ∆x)

     < 0, if x lies inside of T t n+1 I = 0, if x lies on the surface ∂T t n+1 I > 0, if x lies outside of T t n+1 I . (IV.13)
The generalisation of [START_REF] Baerentzen | Signed Distance Computation Using the Angle Weighted Pseudonormal[END_REF] defines three distinct normals that are useful in the signed distance computation, namely the vertex normals, defined using the incident angles α i , as well as the edge normals and the face normals. In fact, the incident angle of an edge is equal to π, so that equation (IV.12) leads to the average

n α xs = πn f 1 + πn f 2 ∥πn f 1 + πn f 2 ∥ = n f 1 + n f 2 ∥n f 1 + n f 2 ∥ . (IV.14)
The incident angle of a face is 2π, so that the equation (IV.12) for the angle weighted pseudo normal simplifies to the face normal itself

n α xs = 2πn f 1 ∥2πn f 1 ∥ = n f 1 ∥n f 1 ∥ . (IV.15)
IV.4 Application of mesh adaptation for an automotive crankshaft The inductor geometry features steep angles and corners. It is modelled using a solid square copper bar that is shaped into the present form. The surfaces are bounded by sharp corners and edges. The bottom part is rounded to enlarge the heat affected zone inside of the crankshaft in the radial direction. This base is directly joined to the elongated connectors resulting in a geometrical description with steep angles.

This complicated inductor geometry can not easily be expressed by analytical functions. Therefore, the discrete triangle surfaces are used to generate a discrete identification function to model the inductor motion. The discrete level set function (IV.9) is based on the correct representation of the pseudo normals, which are computed for the geometric features of each finite element, namely faces, edges and vertices. Each localised computation utilises the incident angle of each feature and can be achieved in parallel.

Figure IV.9 shows the resulting angle weighted pseudo normals for the geometric features of the inductor, as described in equation (IV.12). The pseudo normals are shown in red. It can be seen that the normals on the flat surfaces are pointing in the same directions as the original triangle surface normals. The edge normals, as well as the vertex normals, are formed using the incidence angle weighted normals of the adjacent connecting geometric features. For the top surface it can be noted that the resulting edge normals point in a 45 • angle. The pseudo normals for the vertex point at the corner edges are built using three consolidated surfaces. It follows that these pseudo normals correctly approximate the surface normals of the geometry, so that they can be utilised in the computation of the discrete level set function of the rotated inductor domain.

In order to identify the elements in the enclosing air domain T tn E that are influenced by the movement, the closest feature of each finite element has to be identified. The geometric test has to be done for all geometric features of each finite element. Figure IV.10 depicts an exemplary tetrahedral element, in order to visualise the element identification procedure. The tetrahedron is depicted in blue. As can be seen, its closest feature is the element edge that is nearest to the corner of the rotated discrete inductor domain T t n+1 I

. The angle weighted pseudo normal at the inductor corner can be utilised with the distance vector to the closest feature, shown in red, in order to compute the discrete level set value, given by equation (IV.9). It can be noted that the level set value associated to this finite element is positive, since it is unaffected by the inductor rotation.

Figure IV.11 shows the resulting discrete level set function values for the complete domain for the final rotation at time step t n+1 as a cut in the x-y-plane. The initial inductor domain T tn I is overlayed as a reference. It can be seen that the set of affected finite elements is limited to the area surrounding the initial inductor domain. The negative values of these elements are shown in red. The negative values indicate that these elements are affected by the inductor rotation, i.e. the rotated inductor domain penetrates any geometrical feature of these discretisation elements.

The tetrahedral elements with negative values compose the list of removable elements T

t n+1 E R
that can be safely removed from the enclosing air domain T tn E . The element subset is localised to the inductor movement, so that the finite element mesh can be readapted with standard remeshing schemes. In practice, a simple connecting operation with additional interior vertex seed points can be used to generate the bounded domain. The results demonstrate that the identification step can be restricted to an initial restricting geometry, as proposed in equation (IV.6). This restricting geometry can be used to significantly reduce the computational time, because it reduces the number of possible removable elements based on a simple bounding geometry that reasonably encloses the moving domain.

Possibilities for shapes are detailed in [START_REF] Arvo | Graphics Gems II[END_REF]. These include simple spheres, oriented bounding boxes (OBB) or axis aligned bounding boxes (AABB). OBBs that is related to the principal orientation of the reference geometry. Oriented bounding boxes are simpler to implement, but fixed to the given global coordinate system. Boundary tests that are necessary to identify elements belonging to the subset of possible removable elements are simple to implement. For rectangular shaped geometries all four sides and each of the eight vertex points need to be tested. The computational effort can be minimised by utilising a tree structure that is based on the coordinate directions. Unfortunately, any possible geometric rotation of the underlying geometry can not efficiently be included in the bounding geometry of an OBB, since the global directions of the global coordinate system are fixed. In contrast, AABBs are based on the principal directions of the moveable reference domain itself. Therefore, the bounding geometry can be rotated together with the rigid inductor geometry. The computational effort is slightly elevated, but is outweighed by the reusability of the bounding geometry throughout the full heat treatment computation. shown in orange. It can again be noted that the readapted mesh is highly localised to the surrounding area of the rotated inductor. The element size is evenly distributed, so that the remeshing phase does not introduce distortions. The mesh quality remains comparable to the initial mesh. It can be seen that the localised algorithm results in the adaptation of a very restricted subset of elements, so that the mesh quality remains high.

The local mesh adaptation can now be utilised inside the heat treatment procedure of the complete work piece and inductor system. Figure IV.13 visualises the cut of the final mesh after the initial inductor rotation by 2.5 • . The initial distribution of heat in the crankshaft domain results in an elevated temperature profile close to the crankshaft surface that is nearest to the inductor. The power of the current source is chosen so that the surface temperature of the heat treated zone rises by approximately 20 • C once the next time step is reached. It is visible that the volume of the heat affected area is effectively enlarged in the radial direction, due to the round shape of the inductor bottom part.

The rotated inductor domain at the next time step t n+1 is fully embedded in a precise enclosing domain T t n+1 E , which is coloured in grey. The average finite element size increases closer to the surface area ∂T t n+1 E of the embedding domain. Close to the inductor surface, the embedding elements are smaller and approach the size of the elements of the inductor surface. The details of the fine mesh are clearly visible in the small air strip between inductor domain and crankshaft surface near The orange coloured finite elements indicate the recreated finite element mesh that is generated to fill the void, due to the exclusion of the removable subset. The new elements do not diverge in size or shape from the initial description. The quality of the finite element mesh is comparable to the initial domain and does demonstrably not diminish the discretisation quality, if a reasonably bound on the permissible inductor motion is defined, as demonstrated in equation (IV.3). The highly efficient mesh identification and adaptation algorithm can be used for all subsequent time steps, ensuring the correct coupling between inductor and workpiece.

V Industrial applications

Français:

English:

This chapter details several industrial applications. The goal is to demonstrate the reliability of the numerical model that has been presented in this work. The examples are ordered by increasing difficulty. In the beginning, there is a model of a cylinder that is heated by an inductor in ring-configuration. The results will be compared to the exponential penetration model that has been introduced in a previous chapter.

Afterwards, the numerical model is applied to calculate the induction heat treatment of a gearwheel. It concerns a complex geometry, which necessitates a fine finite element discretisation. The model is treated with the assumption that the material remains constant throughout the simulation. Nevertheless, it is possible to indicate why it is necessary to utilise a numerical simulation to deal with this complex geometry in contrast to an approximation using the exponential decaying model. In effect, the current cancellation effect, which is created in the teeth of the gearwheel, leads to a complex behaviour of the heating power density, which is generated by the eddy currents.

The final example is the automotive crankshaft of project OPTIPRO-INDUX.

It possesses a fully non-linear material behaviour. In the beginning, there is a demonstration of the application of the source current on a complex domain, i.e. the used inductor. Afterwards, the evolution of the temperature is shown, followed by the evolution of the material and a demonstration of the effects that appear. Finally, the good behaviour of the efficient solver is illustrated, which leads to a rapid and stable convergence.

V.1 Cylinder spin hardening

This test case involves a cylindrical shaped part with a diameter of 7.5 mm that is placed in the centre of a ring shaped inductor on which a voltage potential is applied that changes in time using a frequency of 10 kHz. An overview of the model is given in figure V.1. The material properties of the workpiece Ω C are non-linear and depend on the temperature T . The thermal properties are chosen according to an N-alloyed austenitic steel of type 316L (X2CrNiMoN17-12). The electrical resistivity is defined using a certain transition temperature T t , where the character of the behaviour changes. The material behaviour is defined for the temperature T in • C using the following conditional form

if (T < T t ) ρ = ρ 0  1 + α ρ  exp  T τ 1  -1  (V.1) else ρ = ρ 0  1 + α ρ  exp  T t τ 1   1 + τ 2 τ 1  -1  (V.2) ᾱρ = -α ρ ρ 0 τ 2 ρτ 1 exp  T t τ 1 + T t τ 2  (V.3) ρ = ρ  1 + ᾱρ exp  - T τ 2  (V.4) end if,
where α ρ , τ 1 and τ 2 are curvature parameters for the temperature regions before and after the transition temperature T t has been reached and ρ 0 defines a reference electrical resistivity. For the cylinder spin hardening test case, these parameters are chosen as

T t = 200 • C, α ρ = 0.6, τ 1 = 150 • C, τ 2 = 500 • C and ρ 0 = 1.35 • 10 -7 Ωm.
The initial temperature has been defined as T t 0 = 600 • C, so that the material is well above the transition temperature T t . It follows that the electrical resistivity inside the workpiece domain is defined using the non-linear relationship given in equations (V.2)-(V.4). The relative magnetic permeability is defined as

µ r (T ) =  max(µ max r (1 -exp( T -T C αµ )), 1), if T < T C 1, if T ≥ T C . (V.5)
For the cylinder spin hardening test case the material is defined, using a Curie temperature of T C = 850 • C, the curvature parameter α µ = 140 • C and an initial maximum value of the relative magnetic permeability of µ max r = 10. In the following, the evolution of the physical quantities is presented for three distinct time steps. Time step ① is the initial step at the end of the final electromagnetic period, after the convergence of a stable heating power density has been established.

Time step ② is at a later instant of time, after the temperature has been allowed to evolve, but before the Curie temperature has been reached. Afterwards, time step ③ is shown after the Curie temperature has been reached in a large part of the inductor domain.

Figure V.2a gives an overview of the temperature evolution for each time step, measured using a centre cutting line. The initial time step ① gives the magnetic field strength for the final sub-time step of the first heating power density calculation, before the first temperature evolution computation is accomplished. Therefore, the cylindrical part is fully homogeneous, due to the weak coupling procedure of the electromagnetic/temperature evolution problem.

The second time step ② shows an exponentially decaying behaviour, such that the largest temperature is to be found near the surface of the domain, which then decreases in an exponential manner to the centre of the workpiece. It can be noticed that the temperature at the centre of the workpiece is greater than the initial temperature. This is solely due to temperature diffusion effects, since there are no eddy currents in the centre of the material, because of current cancellation effects.

Finally, time step ③ shows a non-exponentially decaying temperature distribution. Its form follows a waveform. This pattern is due to the non-linear interaction of the material evolution and the changes in the applied magnetic field.

The distribution of the relative magnetic permeability is visualised for several distinct time steps in figure V.2b. The relative magnetic permeability decreases with rising temperature. Therefore, the relative magnetic permeability is greatest for the initial time step ①. The second time step shows that the relative magnetic permeability starts to decrease in the complete domain, with smallest values close to the surface of the workpiece. The former time step shows a temperature distribution that is well below the Curie temperature. In contrast, time step ③ shows the relative magnetic permeability for a time step with parts of the domain well in the austenitic region. Therefore, the relative magnetic permeability in large parts of the domain instantaneously jumps to the value of the relative magnetic permeability in free space, i.e. µ r = 1.

The induced heating power density is greatest for the initial time step ①, because the penetration depth is smallest for this time step, as can be seen in figure V.2c. The heating power density is scaled by the maximum value of the first step. The heating power density follows the exponentially decaying form of the approximation that is presented in subsection I.5.1. The induced heating power density of the second time step ② is smaller than that of the initial time step, due to the increase of the penetration depth, as has been mentioned above. Time step ③ shows a highly non-linear behaviour that shows a non-exponential behaviour. In contrast to the approximation using the eddy current distribution of equation (I.6), the maximum value of the heating power is not found near the surface, but in the interior of the domain. Figures V.3a, V.3b and V.3c show the magnetic field for each of the abovementioned time steps with an addition of streamlines, in order to visualise the direction of the generated magnetic field. The cylindrical part, as well as the ringinductor, are shown as solids using the colors grey and red. The skin effect is clearly visible for the first time step in figure V.3a, showing a large concentration of the magnetic field near the surface of the cylindrical part. The magnetic field is symmetric and extends in a ring-like manner in the orthogonal direction of the principal winding direction of the inductor.

The magnetic field strength is largest in this time step, in contrast to time steps ② and ③, as can be seen in figures V.3b and V.3c. This can be explained by the change of the physical properties of the cylindrical part, due to a rising temperature, which results in a decrease of the relative magnetic permeability, leading to an increase of the penetration depth with respect to the applied frequency as has been presented in equation (I.12)

δ ≈ 503  ρ f µ r H -1 2 m 1 2 .
The volume of the affected zone increases, resulting in a decrease of the magnetic field strength. This results in a decrease of the applied heating power density for the later time steps. This effect has been mentioned in [START_REF] Rudnev | Handbook of induction heating[END_REF] (see also [START_REF] Grum | Induction hardening[END_REF]), which gives several examples for the highly non-linear waveshaped pattern of real world induction heating applications. The penetration depth of the cylinder test case is approximately 2 mm, which is larger than the non-magnetic layer visible in figure V.2b. [START_REF] Rudnev | Handbook of induction heating[END_REF] mentions that the wave-like form arrangement is due, in essence, to the remaining magnetic properties inside of the workpiece domain, which is called the dual-properties phenomenon. It is furthermore mentioned that these non-linear effects can lead to an overheating in the interior of the domain, which is not predicted by the approximation theory presented in subsection I.5.1. It is used as a demonstration for the need of a coupled numerical analysis, in order to explain real world heating phenomena. The same wave shaped pattern can be found in [START_REF] Nikanorov | Approaches for Numerical Simulation of High Frequency Tube Welding Process[END_REF], which gives examples regarding high frequency induction welding of steel tubes. Two numerical models are presented. The first model includes a full discretisation of tube, inductor and an impeder. The second model is a localised version including only the welded tube. Both models are defined in a global configuration including a surrounding air domain. The local model is chosen, since it allows to increase the accuracy of the model without exceeding available computing resources. The temperature profile in the welding point of the steel tube (see figure 9 in [START_REF] Nikanorov | Approaches for Numerical Simulation of High Frequency Tube Welding Process[END_REF]) is given for the full wall thickness. It follows a waveshaped curve, which flattens in It should be noted that the complex interaction of thermal and magnetic properties result in a complex system response. [START_REF] Nikanorov | Approaches for Numerical Simulation of High Frequency Tube Welding Process[END_REF] notes that a full model is absolutely necessary for studying the correct system behaviour with respect to parametrical studies and optimisation. Yet, the localised model can be used to study the non-linear effects in the welding point.

It can be concluded that the approximative penetration depth model, using assumptions of exponentially decaying induced eddy currents, is only applicable for simple geometries with fully homogeneous materials.

V.2 Gearwheel spin hardening

This section shows the application of induction heat treatment for the complex geometry of a gearwheel. These workpieces are hardened in the regions of their teeth, mostly on the faces and the flanks. In an industrial setting the heating is usually done in several heating cycles with different frequencies to evenly distribute the heat in the heat affected zone. In the following, the effects of the choice of frequency on the distribution of induced eddy currents is shown for a low frequency of 1 kHz in comparison to a high frequency of 100 kHz. The gearwheel base circle can be identified as the circular shape that would persist if the teeth were cut. This region behaves like a general cylindrical part during induction heating, as described in section V.1. In contrast to the base circle, the teeth show a complicated behaviour, due to current cancellation effects [START_REF] Rudnev | Handbook of induction heating[END_REF]. The magnetic field created by the inductor enters each tooth not only through the top face, but also from both sides of the rest of the flank. This generates opposing eddy currents that cancel each other out, such that the heating power is showing a highly complex, geometry dependent behaviour. In the following, a non-magnetic material for the gearwheel is considered, using a relative permeability µ r = 1 and a temperature independent resistivity ρ = 1.43 • 10 -7 Ωm, such that the current cancellation effects can be visualised without introducing additional effects due to non-linear material behaviour. The thermal properties are chosen according to a carbon steel of type C45.

Figures V.5a and V.5b show the location of two cutting lines that are used to visualise the power density in different subsections of the gearwheel. The wire frame model visualises the finite element model of the workpiece. The first line extends from the centre outwards, cutting only the base circle. The radius of the base circle is 50 mm, with a cylindrically shaped hole of radius 15 mm in the centre.

Figure V.6a and V.6b show the distribution of induced heating energy for an applied frequency of 1 kHz and 100 kHz. The penetration depth for the low frequency test case is δ = 6.015 mm, whereas it is 0.0615 mm for the high frequency test case.

It is visible in figure V.6a that the low frequency test case features a deep distribution of induced eddy currents inside the workpiece. The heating power density following the first line is showing an exponentially decaying behaviour, as predicted by the assumptions on the current distribution using equation (I.6). The penetration depth is indicated and coincides well with the numerical results. The blue-dotted line indicates the heating power density of the second cutting line. It cuts through a tooth. The penetration depth is much larger than the width of each tooth. Therefore, the current cancellation effects lead to a decrease of the induced eddy currents inside the workpiece in these regions. This is clearly visible for the region farther away of the base circle in figure V.6a. The induced heating power density only follows an exponential behaviour inside of the base-circle domain.

Figure V.6b shows the effects of a high frequency on the distribution of induced heating power density, inside the workpiece, using a frequency of 100 kHz. The theoretical penetration depth for this material is indicated in the same manner as for figure V.6a. As predicted by the theoretical assumptions on the eddy current distribution, the heating power density decays rapidly in the direction of the interior domain, such that it practically only heats up the outer diameter of the base circle. For this test case, the penetration depth is smaller than the characteristic diameter of each tooth, such that the current cancellation effects are not very pronounced. The blue-dotted line in figure V.6b clearly shows that the current cancellation effects are weak, such that the heating power density follows an exponential like behaviour, even though it shows some divergences near the transition to the base circle.

For both cases it is visible that the maximum heating power density and, therefore, the maximum induced eddy current is on the outermost point of the base circle.

The distribution of heating power, using the above-mentioned frequencies, leads to an uneven distribution of heat in the workpiece. Even though heat dissipates in the domain, this can lead to an over-heating or under-heating of the teeth. Therefore, the actual heating in an industrial setting has to be done using changing heating cycles including different applied frequencies, in order to evenly apply the heat in the workpiece [START_REF] Rudnev | Handbook of induction heating[END_REF]. Figures V.7a-V.7d show the evolution of the temperature for the low frequency case for different steps of one heating cycle. Figures V.8a-V.8d show the same process with the applied frequency of 100 kHz. For the high frequency test case, shown in figure V.8, it is visible that the heating energy is only applied close to the surface of the gearwheel. The centre of the workpiece and even the centres of the teeth near the roots are mostly heated by conductive effects. [START_REF] Schlesselmann | Numerical Calculation and Comparison of Temperature Profiles and Martensite Mictrostructures in Induction Surface Hardening Process[END_REF] provides several temperature and resulting martensitic phase profiles for the teeth of hardened gearwheels, computed by numerical simulations and derived by experiments. An interesting effect is demonstrated regarding differences in the assumption of the simulation sequence in contrast to the real experiment. [START_REF] Schlesselmann | Numerical Calculation and Comparison of Temperature Profiles and Martensite Mictrostructures in Induction Surface Hardening Process[END_REF] visualises the simulated temperature profile of a gear, indicating a contour shaped hardening pattern. In contrast to the simulation, the experimental result shows a through hardening of the full tooth. [START_REF] Schlesselmann | Numerical Calculation and Comparison of Temperature Profiles and Martensite Mictrostructures in Induction Surface Hardening Process[END_REF] explains that this is due to the neglecting of the necessary holding time between induction heating and quenching procedure. During this production phase, which lasts 0.4 s, temperature diffusion leads to a heat flow in the direction of the core, so that the resulting martensitic structure is different to the simulation.

The results are in accordance with the examples presented in [START_REF] Rudnev | Handbook of induction heating[END_REF]. In addition, [START_REF] Rudnev | Handbook of induction heating[END_REF] shows special flux concentrator designs that can be used to focus the heat effectively on the flanks, without affecting the region of the base circle. These Gap-by-Gap inductors focus the heat either in each unique tooth or in the roots of the gearwheel and are an alternative to the spin hardening using circular shaped ring inductors to create non-uniform contour profiles to enable a wide variety of hardness profiles.

The numerical test case is not fully centred on the symmetrical point of the ring-centre, as can be seen in figures V.7 and V.8. It deviates slightly in the positive x-and y-direction. The impact is larger for the high frequency test case. In essence, this leads to a relative overheating of the teeth in the top region, compared to the lower region. This is a common problem in induction heating, such that homogeneous hardening patterns can only be achieved when the coil is located in full symmetry to the workpiece [START_REF] Rudnev | Handbook of induction heating[END_REF]. Numerical analysis can be used to spot such problems for complex geometrical parts, before the production process begins.

An additional example for the influence of induction positioning is presented in [START_REF] Przyłucki | Induction Heating of Gears -Pulsing Dual-Frequency Concept[END_REF] for the heat treatment of a gearwheel, which mentions that the heating is more uniform if the distance between inductor and workpiece is increased. In order to derive correct temperature fields, [START_REF] Przyłucki | Induction Heating of Gears -Pulsing Dual-Frequency Concept[END_REF] includes a dual frequency heating approach (cf. [START_REF] Grum | Induction hardening[END_REF]). In an initial phase the workpiece is heated by applying a medium frequency of 10 kHz for 1.5 s, followed by a high frequency heating phase using a frequency of 100 kHz for 0.2 s.

An alternative approach consists in the application of superimposed frequencies (see, e.g. [START_REF] Biasutti | Induction hardening of complex geometry and geared parts[END_REF], [START_REF] Kobos | Analysis of Dual-Frequency Generators Structures for Inductor Heating[END_REF] or [START_REF] Schwenk | Experimentally validated residual stresses and distortion prediction for dual frequency induction hardening[END_REF]). This approach is based on a direct application of several superimposed current sources. The simultaneous application is motivated by noting that it can then be avoided to switch between frequencies or distinct coils during the heat treatment [START_REF] Biasutti | Induction hardening of complex geometry and geared parts[END_REF]. The heating cycle is shorter, so that the possibility of cool-downs or unwanted temperature diffusion is minimised. [START_REF] Schwenk | Experimentally validated residual stresses and distortion prediction for dual frequency induction hardening[END_REF] gives an experimental validation of a dual frequency induction heating approach, with a focusses on the residual stress profiles in the hardened zone, due to the material phase transformation, concluding that induction heating using superimposed middle and high frequency sources is a promising approach.

[ [START_REF] Biasutti | Induction hardening of complex geometry and geared parts[END_REF] presents several examples with respect to the dual frequency heating of a gearwheel. It notes several phenomena, which affect the distribution of the induced eddy currents. It mentions that the proximity effect results in nonuniform distributions inside the inductor. Furthermore, it is mentioned that the end effect and edge effect result in non-uniform magnetic fields in the axial direction of the heat treated workpiece, resulting in a variation of the induced heating power in the axial direction. As a remedy, it is found that the coil size can be adapted to counteract this effect and that flux concentrators should be placed on top and below the workpiece. The edge effect results in an increase of the induced heating power in the top region of the teeth near the corners, due to an interaction of the axial and radial components of the magnetic field at these points. [START_REF] Biasutti | Induction hardening of complex geometry and geared parts[END_REF] mentions that it is affected by the applied heating frequency, as well as the coil geometries.

Additional examples for the influence of the edge effect on the temperature evolution of heat treated gear wheels are presented in [START_REF] Bocher | Simulation of Fast Induction Surface Heating and Comparison with Experimental Full-Field Surface Temperature Measurements[END_REF] (cf. [START_REF] Candeo | Multiphysics Modeling of Induction Hardening of Ring Gears for the Aerospace Industry[END_REF]). The experimental setup involves a high frequency power supply provided by EFD Induction, which generates source currents with a frequency of 190 kHz. [START_REF] Bocher | Simulation of Fast Induction Surface Heating and Comparison with Experimental Full-Field Surface Temperature Measurements[END_REF] shows comparable results to [START_REF] Biasutti | Induction hardening of complex geometry and geared parts[END_REF], demonstrating a relative overheating of the region of the teeth. It mentions that the edge effect is, in essence, influenced by the electromagnetic material properties of the model, whereas the temperature gradient is more influenced by the thermal parameters. It is shown that many physical effects can be reproduced, but it concludes that the correct reproduction of the full thermal history of a surface heating procedure is still an open problem.

[ [START_REF] Kobos | Analysis of Dual-Frequency Generators Structures for Inductor Heating[END_REF] analyses dual frequency generators for induction heating applications. It gives several examples of output waveforms generated by different generator sources (see, e.g. figure 2 and figure 3 in [START_REF] Kobos | Analysis of Dual-Frequency Generators Structures for Inductor Heating[END_REF] or figure 8 in [START_REF] Biasutti | Induction hardening of complex geometry and geared parts[END_REF]).

The numerical model proposed in this work is based on a fully transient formulation of electromagnetic vector potential. The applied frequency can, therefore, be applied using an arbitrary shape, which would make it particularly useful for the above-mentioned dual frequency simulations using superimposed source currents.

The finite element mesh of this complex geometrical part includes around 1.4 million elements, resulting in around 4 million edge degrees of freedom for the electromagnetic problem. The linear system of equations for the electromagnetic problem can be represented using slightly over 65 million non-zero entries. The convergence is uniform throughout the computation, as can be seen in figures V.9a and V.9b, even though the geometry is highly detailed, featuring steep edges, many possible re-entrant corners and sharp changes between element sizes. The solution technique includes a GMRES solver using an auxiliary space one-cycle algebraic multigrid preconditioner, including a CLJP coarsening technique with a coarsening strength threshold of 0.25. This configuration is chosen, because the material is defined as linear and homogeneous, such that the initial electromagnetic computation can be used throughout the full heating cycle. The blue line in each respective figure shows the convergence of the fastest converging time step, whereas the red line shows the time step that needed the most iterations to reach convergence. The residual error has been scaled to the same initial value, in order to enable a direct comparison of the convergence behaviour.

The convergence is mostly uniform and depends strongly on the frequency of the applied loading. This is due to the fact that the frequency changes the admissible time step size, i.e. it decreases for higher frequencies. Therefore, the linear system of equations, as well as the right-hand side, are weighted differently with respect to the low frequency test case.

Figure V.9a shows that convergence can be achieved in a very low number of iterations for the low frequency test case. On average, the linear system of equations converges in 24 iterations.

The high frequency test case converges weaker, needing on average 51 solver iterations to reach convergence. Yet, both cases show that the convergence is very uniform with only slight variations for the different loading cycles, since both the median and the mean are equivalent for each case. An overview of the computational time for the electromagnetic computation is given in table V.1 for the 1 kHz, as well as for the 100 kHz configuration using an auxiliary space multigrid preconditioned GMRES solver in connection with a CLJP coarsening technique in a sequential configuration. This is in accordance with the numerical experiments by [START_REF] Sterck | Reducing Complexity in Parallel Algebraic Multigrid Preconditioners[END_REF], which shows that the HMIS convergence technique converges much slower than the CLJP technique. Furthermore, it shows that the convergence behaviour of the HMIS technique is severely affected by the granularity of the coarsening. For a 3D elliptic problem with varying material coefficients it is shown that the HMIS technique needs approximately 12 times more solver iterations to converge than a CLJP technique for the same granularity for a test case with a distribution of approximately 64.000 nodes per processor. For the gearwheel test case, it can be observed that the strong loss of convergence is observed for the 16 processor test case, which results in a rough distribution of 200.000 degrees of freedom per processor. The setup time for each coarsening technique is visible in figures V.10b and V.10e. It can be seen that the setup time peaks for the 2 processor test cases, since the effects of the parallel computation, i.e. parallel messaging due to process communication, together with a large coarsening surface, lead to a large increase in computational time. The curve decreases for a rising processor count and flattens for the highly parallel test cases. The good convergence behaviour of the CLJP technique leads to lower solution times, as can be seen in figure V.10c. It can be noted that the efficiency is greatest for the 8 processor configuration. An increase of the granularity leads to an increase in process communication and a general loss in numerical efficiency. This can also be seen in figure V.10f for the HMIS coarsening technique. For this case, it is visible that the bad convergence behaviour leads to a strong increase in the solution time. In comparison to the CLJP technique, the most efficient parallel configuration is achieved by using a 4 processor setting. The applied electrical loading varies with a frequency of 10 kHz. The inductor is modelled using a fully closed domain, in order to achieve conformity with the weak form of the electromagnetic problem, such that uniqueness and fast convergence is ensured. The voltage potential distribution on the inductor domain is visible in figure V.12. It can be seen that the inductor is modelled with sharp edges, especially at the points where the generated electrical current changes directions. This is due to using a simplified CAD geometry that only approximates the form of the real inductor. Sometimes, these design inaccuracies can lead to convergence problems, which makes this inductor a reasonably problematic boundary test case. It will be shown later in this section that the convergence is, in fact, not negatively influenced by these sharp edges.

The flux concentrators, visible in black, are placed on top and inside the loop, created by the inductor wiring. They are used to focus the magnetic fields. Possible materials are ferrites or other metals or ceramics with a high relative magnetic permeability. The flux concentration effect is shown later in this section. The initial computation of the voltage potential on this complex domain is achieved using the double cutting plane technique that is shown in subsection II.3. Figure V.13 shows the 4 distinct cutting surfaces for this inductor. The cutting plane on top with normal in the x-direction results in two separate surfaces marked in light-blue and orange. The second cutting plane in the y-direction, which is also visible in figure V.12 creates two additional cutting surfaces, shown in blue and green.

The initial cutting plane is used to compute an initial voltage potential distribution using the formulation (II.58)-(II.60). The elements are not aligned perfectly in the direction of the cutting plane, thus the strong imposition of the Dirichlet conditions leads to slight inaccuracies near the cutting plane. The current at this section will be slightly divergent. The second cutting plane is chosen in a straight section, where it can be assumed that the computed current of the first computation will be uniformly pointing in the y-direction. The second cutting plane can therefore be used with the resulting nodal values of the initial computation, in order to compute a fully conforming voltage potential, resulting in a smooth non-diverging source current.

Afterwards, the source current is projected using the curl-conforming potential representation presented in subsection II.3, in order to ensure the full conformity of the applied source current. This procedure can be used for arbitrarily complex inductors, to ensure the conformity of the applied currents and the fast convergence of the electromagnetic magnetic vector potential problem. The electromagnetic time step shown is the final step of the first electromagnetic computation with a homogeneous relative magnetic permeability, before the convergence of the heating source density is established. Therefore, this test case is equal to a fully linear homogeneous test case.

It can be seen that the induction application without flux concentrators, i.e. a relative magnetic permeability equal to the one of the surrounding air-domain, leads to a wide distribution of the magnetic field in the crankshaft domain. Figure V.14c shows the homogeneous distribution over the circular domain.

In contrast, figures V.14b and V.14d show the effect of using flux concentrators with a high relative magnetic permeability on the computed magnetic field B. The magnetic field is concentrated below the flux concentrators. In particular, this effect is visible in figure V.14d, which shows that the flux concentrators lead do a large relative decrease of magnetic field strength near the centre. The concentration of the magnetic field throughout the computation is greatest for the case using a flux concentrator with a relative magnetic permeability µ r = 1000, even though the concentration effects of the test case with µ r = 100 is only slightly smaller. The main difference is visible in contrast to the case shown in figure V.15a that results in a very homogeneous distribution of heat under the inductor. The application volume of the heating power is enlarged, such that the final maximal temperature is much lower than for the cases involving efficient flux concentrators.

In the following, the fully non-linear effects of the numerical model will be demonstrated for the test case using flux concentrators with a relative magnetic permeability of µ r = 100. The material model of the workpiece shows a temperature dependent behaviour. Its thermal properties are defined according to a steel of type 38MnSi4. The relative magnetic permeability is described using the model presented in equation (V.5), using a Curie temperature of T C = 732 • C, the curvature parameter α µ = 80 • C and an initial maximum value of the relative magnetic permeability of µ max r = 10. The electrical resistivity is modelled using equations (V.2)-(V.4) with given transition temperature T t = 200 • C. The reference resistivity is given as ρ 0 = 1.35 • 10 -7 Ωm, whereas the curvature parameters are defined as α ρ = 0.6, τ 1 = 150 • C and τ 2 = 500 • C. The non-linear material models for the electromagnetic parameters are visualised for the given temperature range of this test case in figures V.16a and V.16b. Figure V.17a gives an overview of the temperature distribution after the material has been allowed to evolve. The temperature range is slightly above the Curie temperature in some parts of the domain. It can be seen that the temperature evolution is not only influenced by the position of the flux concentrators, but also by the geometrical features of the crankshaft itself; like the hole used for fixing the bearings.

The relative magnetic permeability for this time step is visualised in figure V.17b. The relative magnetic permeability is mostly homogeneous in the crankshaft domain. It decreases only in the region close to the inductor, where the temperature is elevated. The blue regions indicate the domains where the Curie temperature T C is reached, such that the material is in austenitic configuration.

The electrical conductivity σ, seen in figure V.17c, closely resembles the distribution of the relative magnetic permeability. Yet, it can be seen that its distribution is smoother, due to the fact that there are no discontinuities in the material evolution curve V.16a.

In the following, four different time steps are shown for a cut in the x-y-plane directly under the flux concentrator closest to the origin. Figure V.18 shows the cut for the first final converged electromagnetic computation, defined as ①. The homogeneous temperature is shown in V.18a, whereas V.18b shows the relative magnetic permeability. The line, indicated in V.18b is used to enable a comparison of the temperature, the relative magnetic permeability, as well as the heating power density in each time step. Figures V.19a and V.19b show the temperature and relative magnetic permeability in a later time step, indicated as ②. The temperature has increased, but remains well below the Curie temperature, such that the relative magnetic permeability, shown in figure V.19b, is well above µ r = 1 in the whole domain. Figures V.21a and V.21b show the temperature and relative magnetic permeability for the time step ④, well after the Curie temperature has been reached. The large domain, shaded in dark-blue in figure V.21b is an indication that large parts of the crankshaft domain under the flux concentrator have been transformed into the austenitic phase. The temperature distribution in figure V.21a shows that the temperature on the sides has increased well above T C , whereas the centre remains relatively cold. This is due to the fact that the induced heating power density is zero at the centre, because of current cancellation effects. The temperature can only increase through conductive effects.

The overall temperature range begins from T = 593 • C, since it is based on the temperature in the full crankshaft domain. Radiation effects and the exchange between the crankshaft surface and the surrounding air, have let to a cooling in parts of the domain that have not been influenced by the induced eddy currents.

The cutting line, visible in figure V.18b enables a comparison of the heating power density for each time step. The comparison is shown in figure V.22. Each heating power density is scaled by the maximum value of the first time step ①. The heating power density for the homogeneous temperature and, therefore, the homogeneous distribution of relative magnetic permeability µ r is shown as a black .19a. The penetration depth is slightly increased, such that the maximum value of the heating power density is lower than for the fully homogeneous test case.

The red line indicates the heating power density for time step ③, shortly after the Curie temperature has been reached in parts of the domain. The penetration depth has increased further, resulting in a decrease of the maximum value of the heating power density and a deeper distribution of heat in the crankshaft domain.

Finally, time step ④ is indicated in yellow for a point in time, where a large domain of the heat affected zone has reached the Curie temperature. The penetration depth is largely increased, such that the heat is applied deeper in the material than for the former time steps. Yet, the large increase in penetration depth results in a large reduction of the maximum value of the applied heating power density. The finite element model, including the global surrounding air domain, is discretised using around 2.6 million elements, resulting in a linear system of equations for the electromagnetic vector potential computation, consisting of around 3 million edge degrees of freedom. It is described using roughly 50 million non-zero entries.

The linear system of equations is solved using a preconditioned GMRES solver, with an auxiliary Maxwell space one-cycle algebraic multigrid preconditioner, featuring an HMIS coarsening technique and a coarsening strength threshold of 0.25. In contrast to the CLJP coarsening scheme, the HMIS coarsening technique converges slower, but is easier to set up, as has been shown in chapter III. Therefore, the HMIS coarsening technique is chosen, because the non-linear behaviour of the material necessitates a reassembling of the linear system of equations during the temperature evolution. The disadvantage is the weaker convergence of this coarsening scheme, in relation to the stronger converging methods.

Figures V.23a, V.23b and V.23c show the convergence behaviour for a full numerical solution for three different test cases, using no flux concentrators in the test case shown in figure V.23a, using flux concentrators with a relative magnetic permeability of µ r = 100 in figure V.23b and using flux concentrators with a relative magnetic permeability of µ r = 1000 in figure V.23c. The relative residual error is scaled to the same value, in order to enable a direct comparison.

Each figure shows the fastest converging time step for the full non-linear analysis in blue, whereas the slowest converging time step is indicated in red. It is visible that the convergence rate is very homogeneous, since the average value is very close to the median in each test case. Convergence for each case can be achieved in slightly over 50 solver iterations.

Interestingly, the large change of relative magnetic permeability µ r has no profound impact on the convergence behaviour, i.e. the solution time is not highly influenced by the choice of relative magnetic permeability of the flux concentrators. It is visible that the test case using a flux concentrator with µ r = 100, visible in figure V.23b, shows a slightly worse convergence than the case with no flux concentrators, visible in figure V.23a. Yet, the test case with the flux concentrators featuring the highest permeability of µ r = 1000, visible in figure V.23c, converges faster than both former test cases.

It is difficult to assess the reason behind this phenomenon, since the coupled analysis features distinct temperature dependent material models that each influence the distribution of the magnetic field, due to the impact on the penetration depth. Judging from the results of figures V.23a-V.23c it is, therefore, reasonable to assume that the choice of relative magnetic permeability for the flux concentrators has no profoundly negative impact on the computational time.

The non-linear material evolution has been demonstrated, as well as the effect of changing penetration depths on the distribution of the heating power density. The computational representation of the crankshaft involves many geometrical features that can be considered to be problematic with respect to numerical computations, like holes, steep edges, sharp corners, larger material coefficient jumps and large differences in relative finite element mesh sizes.

This test case shows the applicability of the preconditioned GMRES solver using the auxiliary space algebraic multigrid method, for solving the global linear system of equations of the electromagnetic vector potential formulation. The introduction of conforming source currents for the complex inductor geometry, using the dual cutting plane approach, in connection with the curl-space projection of the source current density, results in fast convergence with a large reduction of relative residual errors. This test case shows the feasibility of the weak coupling approach for the induction heating simulation for complex workpieces.

An overview of the fully coupled electromagnetic/temperature evolution computation is given in table V The good convergence behaviour of the auxiliary space Maxwell algebraic multigrid solver in connection with a CLJP coarsening technique is visible in figure V.24a. For this test case the convergence behaviour is not affected by the parallel coarsening of this method. In contrast, figure V.24d shows the convergence behaviour using a HMIS coarsening. It is visible that the convergence behaviour is better than for the geometrically complex gearwheel test case, shown in figure V.10d. It can be assumed that this effect is due to the more homogeneous distribution of elements and materials in the crankshaft test case. The gearwheel test case features strong changes in element size, coupled with the highly detailed discretisation of the teeth. The coarsening leads to a dispersion of the initial gearwheel surface that is increased with rising granularity. In contrast, the crankshaft test case features small elements in the region directly between inductor and crankshaft surface. The element size decreases smoothly towards the outside surface of the air domain, thus leading to a more equal distribution in the coarsening step. Even though the crankshaft test case shows a better convergence behaviour than the gearwheel test case, the negative effects of the parallel partitioning are also visible. Interestingly, the highly parallel test case using 32 processors converges slightly better than the test case using 24 processors, as can be seen in figure V.24d. This is an indication of the complex relationship between coarsening, parallel partitioning and numerical convergence.

The parallel setup time, visible in V.24b and V.24e follows the same behaviour as the gearwheel test case V.10b and V.10e. Again, the setup time is related to the surface coarsening area and the amount of parallel coarsening. The setup time peaks for the 2 processor configurations, decreases with increasing parallelism and flattens out for higher processor counts. The advantage of the HMIS coarsening technique should rest in the faster setup times. Interestingly, this effect can only be noted for the sequential, as well as for the highly parallel test cases. Unfortunately, the gains in the highly parallel configurations are marginalised by the loss of efficiency in the solution step.

The solution time shows the influence of the convergence, as well as parallel effects, like reduced efficiency due to increased inter-processor communication. The CLJP solver is more efficient, due to its better convergence. The disadvantage of the HMIS coarsening technique is clearly visible for the highly parallel test cases, shown in figure V.24f. The increase of the processor count leads to an increase in computational time. The same effect is visible for the CLJP coarsening technique, even though it is not as strong as for the HMIS coarsening. For the HMIS coarsening it is shown that for this configuration a 4 processor partitioning leads to a faster solution than both a 24 and a 32 processor partitioning. Figure V.24c shows that the fastest solution can be attained using a CLJP coarsening using a 16 processor partitioning. The 24 processor test case necessitates approximately the same computational time as an 8 processor test case. V.3 shows the reduction of the residual error with respect to the solution time for the first increment of the crankshaft test case in a parallel setting using 16 processors. The computational time increases linearly, with the reduction of the residual error. It can be seen that a reduction of the numerical error by 10 -10 is only twice as costly as the reduction of the residual by 10 -5 . The results indicate that the auxiliary space algebraic multigrid technique is useful in applications, which have a need for a large reduction of the residual error.

residual error reduction 10 -5 10 -6 10 -7 10 -8 10 -9 10 -10 time in s 29.

VI Conclusions

Français:

Le dernier chapitre résume les principaux résultats. Les hypothèses concernant le modèle numérique sont comparées avec les conclusions sur les applications de référence, afin d'évaluer les domaines d'application possibles des méthodes numérique proposées. Les résultats acquis concernant les méthodes de résolution efficaces sont présentés pour le traitement numérique des applications industrielles à grande échelle. Une conclusion est ensuite donnée au sujet de l'algorithme proposé pour l'adaptation de maillage, y compris certains domaines d'application supplémentaires possibles. Enfin, les perspectives pour la recherche future et des suggestions pour l'amélioration du modèle numérique sont données.

English:

The final chapter sums up the key results. The assumptions regarding the numerical model are compared to the findings of the benchmark applications, in order to assess possible scopes of application of the proposed numerical method. The findings regarding efficient solution methods are presented in view of the numerical treatment of large scale industrial applications. Afterwards, a conclusion is given regarding the proposed mesh adaptation algorithm, including some possible additional application areas. Finally, an outlook for future research and suggestions for the improvement of the numerical model are given.

VI.1 Numerical model and industrial applications

The numerical model for induction heating involves electromagnetics, heat diffusion and material phase changes. Each of these sub-problems interacts within a global framework that needs to be efficiently coupled. In this context, the electromagnetic formulation has been identified as the most difficult problem from a numerical perspective. The nature of induction heating allows to model the problem using a partial differential equation of parabolic type, due to the quasi-steady state of the electromagnetic induction process. The resulting induced eddy currents can be post-processed to a heating energy term that can be included in a coupled heat diffusion formulation. The parabolic nature of the heat diffusion equation ensures that the coupled system is solvable using standard iterative solution techniques, leading to smooth temperature distributions.

The ultra-weak coupling procedure has been identified as advantageous, since it allows to introduce a coupling between thermal effects and electromagnetic computation involving different time scales. The bridging of the different time scales allows to focus on the electromagnetic problem and the heat diffusion problem separately. The electromagnetic/thermal coupling procedure is based on an averaging of the induced eddy currents on the electromagnetic time scale. The computation is controlled, by user supplied parameters that allow to influence the strength of the coupling, allowing a fine grained control of the approximation of the non-linear material evolution. The parameters define thresholds for admissible material changes. The weak coupling results in a re-evaluation of the thermal source terms, after the relative change of a material quantity reaches the prescribed thresholds, ensuring a correct approximation of the heating source contribution.

The numerical model allows for spatially varying material coefficients, with the possibility to include a vanishing electrical conductivity σ in parts of the domain. Both the electrical conductivity σ, as well as the relative magnetic permeability µ r are allowed to vary with respect to both the temperature and the magnetic field. An algorithm for the computation of source currents on arbitrary closed inductor domains is presented. It ensures the correct introduction of the source currents, so that they are divergence free, resulting in conforming source current distributions leading to a fast reduction of the residual error during the solution stage. The fully transient application of the source current distribution on arbitrary complex closed inductor geometries allows the application of non-linearities on the electromagnetic time-scale. The source currents need not to be sinusoidal so that it is feasible to include more complicated excitation currents. This opens up the possibility to directly model the full induction circuit system including the system response of the power source (see, e.g. [START_REF] Blinov | Mutual influence between load and frequency converter in the induction heating of steel[END_REF] or [START_REF] Dughiero | Coupled Multiphysics Circuital Modelling of Quasi-Resonant Induction Cooktops[END_REF]).

Several example applications have been presented, in order to show the validity of the approach. Large scale industrial examples are provided, including the surface heat treatment of an automotive crankshaft, as well as the spin hardening procedure of a gearwheel. The applications show the viability of the numerical model in connection with efficient solution methods. The automotive crankshaft test case demonstrates the impact of non-linearities during the solution phase, resulting in complex heating power and temperature distributions. The non-linear behaviour of the electrical conductivity and the relative magnetic permeability lead to chan-ging penetration depths. These non-linearities affect the temperature distribution, which in turn affects the generated magnetic field. The results are compared to the theoretical skin depth formula, validating the numerical model.

The gearwheel spin hardening test case shows the advantages of using a numerical model for computing the temperature distribution for a problem involving a relatively complex geometry. The impact of different frequencies on the skin depth is shown. Even though the material properties are homogeneous and chosen to be constant, with respect to the applied magnetic field and temperature. It can be seen that the current cancellation effect leads to a complicated distribution of the induced eddy currents, so that it does not suffice to use simpler methods like the skin depth approximation formula.

VI.2 Efficient linear solvers for the associated electromagnetic problem

The electromagnetic magnetic vector potential problem is implicitly discretised in time, allowing for arbitrary source currents. The applications presented in this work introduce a source current with sinusoidal excitation that is discretised in time using equally spaced time steps in connection with an implicit Euler discretisation.

The linear system of equations is sparse, symmetric and semi definite with a large null space. It is shown that the classical Krylov subspace iteration techniques can be employed for the repeated solution of this linear system of equations. It is demonstrated that the GMRES method is superior to the CG, CGS, BICG and BICGStab methods. The advantage of this solver type is its monotonic convergence behaviour, so that the residual error can be steadily reduced with increasing number of iterations. The GMRES method can be effectively employed if the preconditioner is efficient, so that it is not necessary to restart the method by discarding the generated Krylov subspace.

It is shown that the classical preconditioners, like the Jacobian, the SSOR or the IC preconditioning techniques are relatively ineffective. In that regard, it is demonstrated that the correct treatment of the large null space is the key to generating efficient solution techniques. The auxiliary space preconditioning technique is shown to lead to very effective preconditioning techniques. It is based on the description of the solution vector with respect to a nodal auxiliary finite element mesh. The residual error of the solution vector that is discretised using an edge element approach can very effectively be reduced, after it has been projected onto this nodal auxiliary space with subsequent application of complying preconditioning operators. The transformation operators can be generated with the existing adjacency information of the edge finite element discretisation. The bilinear forms of the resulting finite element decomposition show an elliptic behaviour. This allows the efficient application of standard multigrid techniques.

Since the auxiliary space is based on a nodal description, it is very advantageous to apply an algebraic multigrid coarsening procedure to generate efficient preconditioners. It is only necessary to employ approximative solution techniques. The combination of the GMRES method in connection with a one-cycle algebraic mul-tigrid preconditioner using an auxiliary space decompositioning technique leads to very efficient solution schemes that enable the computation of large scale real world applications. An important aspect in that regard is the choice of the coarsening technique for the discrete operator. To that extent several up-to-date coarsening techniques are presented, including the RS3, CLJP, Falgout, PMIS, HMIS and the ECGC technique. The impact of each coarsening method is demonstrated for large scale industrial applications in sequential and parallel configurations.

A detailed investigation regarding the effects of different solver configurations is presented for an automotive crankshaft test case. The source current is applied using a closed inductor of complex geometry. The finite element mesh includes around 3 million finite elements, resulting in around 3.5 million edge degrees of freedom. The sparse linear system of equations consists of approximately 55 million non-zero values. Additional engineering applications have been demonstrated for a detailed model of an industrial crankshaft with approximately 3 million edge degrees of freedom, resulting in around 50 million non-zero entries and for a gearwheel test case with approximately 4 million edge degrees of freedom, resulting in around 65 million non-zero entries.

It is shown how the coarsening threshold affects the set up, as well as the solution times, since it effects the operator complexity of the coarsened operator. In that regard it is demonstrated that a high approximation quality for the operator related to the large kernel of the linear system of equations is indispensable for a high error reduction rate.

It is found that the CLJP and Falgout techniques lead to coarsened operators that can very effectively be used in sequential or lowly parallel configurations. A large operator complexity leads to a low number of iterations, during the solution stage.

For highly parallel test cases and problems with changing material coefficients, which might necessitate a recreation of the preconditioning operator, it is shown that the HMIS, as well as the ECGC method can be employed very efficiently. In comparison to the CLJP, Falgout and especially the RS3 method, the setup times are reduced. The lower complexities lead to an increase in the iteration numbers, due to a decrease of the convergence rate, but the lower cost of application leads to an efficient scheme in a highly parallel setting. Yet, it must be noted that test cases are provided for which an increase in iterations leads to a large increase in computational time, which offsets the lower application cost.

It is furthermore shown that changing material properties, including many discontinuous effects, do not lead to a large increase in the solution time. This shows the versatility and robustness of the presented solution scheme.

The solution of the linear system of equations can be handled in parallel, so that scalability and an efficient solution time is ensured. Yet, it must be noted that a simple parallel application might not directly translate into a reduction of computational time. In several example problems it is shown that the parallel efficiency is reduced, if the local grid size becomes small. The result is an increase in solution time, which can be even higher for highly parallel configurations, compared to lowly parallel configurations. This is a direct result of inter-processor communication and certain difficulties related to the correct treatment of the coarsening operators for nodes on processor boundaries. In a parallel setting, the initial grid is divided into several distinct parts. An application of the parallelised operator necessitates parallel communication over processor boundaries. An increase in parallel granularity, i.e. an increase in the number of processors, results in the decrease of the average local grid size and, therefore, leads to a growth of the inter-processor boundaries. In addition, it must be noted that the auxiliary space preconditioning technique is based on a nodal auxiliary subspace. It follows that the bandwidth of the algebraic multigrid operator should be optimised with respect to the nodal adjacency numbering. Yet, the matrix bandwidth in this work is only optimised for an edge element approach, so that the discrete elliptic operators of the associated bilinear form of the preconditioning operator might possess a large bandwidth, which could lead to an unnecessary increase in inter-processor communication.

VI.3 Modelling inductor motion

Another important aspect with respect to the numerical treatment of industrial induction heating problems is the introduction of a relative movement of workpiece and inductor domains. A novel method has been presented, which is based on the rigid approximation of the moving inductor domain. The resulting algorithm is based on the discrete description of this movement, using a level set function that can be generated using the surface primitives of the rigid inductor domain. The novel algorithm is very efficient and it is based on local computations, so that it can effectively be computed in a parallel configuration. The algorithm depends on the identification of the relative positioning of discrete geometrical elements like points, lines and faces with respect to the discrete inductor geometry, which is defined using surface elements. The angle weighted pseudo normal has been used to correctly identify this relative position in a very efficient manner, which is valid for inductors of arbitrary shape. It can be noted that the resulting algorithm is independent of the element size, element type and polynomial order, so that it can be used for problems involving tetrahedral or hexahedral formulations of any shape. The application is controlled by a user supplied constant that bounds the maximum relative change of the inductor positioning. The transformation can be monitored during the computation. If it exceeds the prescribed permissible range, the time step size can be reduced until the transformation results in a conforming displacement.

An application for a heat treated automotive crankshaft including the rotation of a complex inductor is provided. The discrete level set identification technique is utilised to readapt the global finite element mesh. The recreated mesh is presented for a rotation of 2.5 • . The identified removable subset of elements is minimal and restricted to the embedding elements closest to the rotating inductor. It is shown that the readapted mesh possesses a very high quality that is comparable to the initial finite element mesh.

VI.4 Outlook

It is shown that the numerical solution techniques involving the auxiliary subspace preconditioning technique lead to a large reduction of computational time. This opens up the possibility to compute large scale numerical simulations, involving highly detailed numerical models of real world problems. To that extent, the numerical model needs to be validated by industrial heat treatment test cases.

The applications in this work involve material non-linearities with respect to changing temperature fields. The influence of the magnetic field on the relative magnetic permeability has not been included in these test cases. This can result in approximation errors with respect to the heating source contribution, due to an inaccurate reproduction of the correct penetration depth of the induced eddy currents. Nevertheless, the numerical model allows the introduction of a dependency with respect to the changing magnetic field following from the computation of the electromagnetic magnetic vector potential. It should be quantified how large the effect of the non-linear material behaviour is with respect to an industrially sized workpiece, like an automotive crankshaft. Additional effects like the holding time between induction and quenching procedure must be included in the simulation, since it has a profound effect on the resulting martensitic phase profiles, as has been demonstrated by [START_REF] Schlesselmann | Numerical Calculation and Comparison of Temperature Profiles and Martensite Mictrostructures in Induction Surface Hardening Process[END_REF].

In addition, the fully transient time discretisation is achieved using an implicit Euler time discretisation, which might be replaced by higher order singly diagonally implicit Runge Kutta (SDIRK) methods to increase the numerical accuracy.

The source currents in this work are modelled using sinusoidal forms. More complex source excitations are possible, which include the direct system response of the induction power supply, e.g. it could be advantageous to include the full transistor circuits, which might generate source currents containing harmonics or to introduce currents with non-sinusoidal waveforms. This would allow a greater control over the heat affected zones and might lead to smoother temperature profiles, as shown by [START_REF] Biasutti | Induction hardening of complex geometry and geared parts[END_REF]. An example for a coupled finite element model between a quasi resonant (QR) converter power source and an electromagnetic induction heating model is given in [START_REF] Dughiero | Coupled Multiphysics Circuital Modelling of Quasi-Resonant Induction Cooktops[END_REF]. This coupled model is an example for a fully transient analysis, including non-sinusoidal waveforms (see figure 7 in [START_REF] Dughiero | Coupled Multiphysics Circuital Modelling of Quasi-Resonant Induction Cooktops[END_REF]). It is used to model the real system response. It mentions that a classical ICCG solver results in unreasonable simulation times, which is why a parallel direct solver has been utilised to reduce the simulation time by a factor of approximately 14. Alternatively, the auxiliary space multigrid preconditioners could be used to increase the scalability of the numerical model.

The proposed mesh adaptation algorithm, which is based on a discrete level set function, is applicable to different problems, for which the movement can be described in a Lagrangian setting. The precondition is that the simulation includes a moveable rigid geometrical feature. Possible example applications are adaptations to existing Lagrangian formulations which necessitate changes in the model geometry, like the replacement or change of rigid walls or the displacement of rigid tools in a forming or forging simulation.

Another important aspect for induction heating design is process optimisation. Many simulation parameters can influence the outcome of a simulation. E.g. the shape, relative placement of the inductor or the frequency and magnitude of the applied source current can be changed. The optimisation process can be reformulated as the minimisation of an objective cost function with different constraints (see, e.g. [START_REF] Favennec | The ultraweak time coupling in nonlinear multiphysics modelling and related optimization problems[END_REF], [START_REF] Rapoport | Optimal Control of Induction Heating Process[END_REF], [START_REF] Di Barbara | Multi-Objective Optimization for Advanced Design of Induction Through Heaters[END_REF] or [START_REF] Touzani | Mathematical Models for Eddy Currents and Magnetostatics[END_REF]). Typical objectives in an industrial application are the optimisation of the energy consumption, the control of the final temperature or to ensure that the temperature does not exceed a given value during the treatment. The optimisation procedure entails a search in the space of possible solutions, so that many sample points need to be computed, resulting in large requirements for computing power. It has been demonstrated that the fully transient magnetic vector potential formulation can be used in the context of a weak coupling procedure to compute industrially sized problems in a reasonable amount of time. The efficient treatment of the resulting semi definite linear system of equations is ensured by the application of the auxiliary space multigrid method. In the future, the model might be applied to optimal control problems, such as the optimisation of the heat affected zone of the presented automotive crankshaft. ABSTRACT : Electromagnetic induction heating is an efficient process allowing to directly heat up a prescribed area beneath the surface of metallic workpieces to enable quenching. This work presents a mathematical model for the coupled electromagnetic/heat transfer process as well as numerical solution methods. The electromagnetic model is based on a magnetic vector potential formulation. The source currents are prescribed using a voltage potential formulation enabling the modelling of arbitrary inductor geometries. The heat transfer problem is modelled using the classical heat diffusion equation. The electromagnetic model is fully transient, in order to allow the introduction of non-linear effects. The space discretisation is based on an edge finite element approach using a global domain including air, workpiece and inductor. The resulting linear system of equations of the implicit formulation is sparse and semi-definite, including a large kernel. It is demonstrated that a preconditioner based on the auxiliary space algebraic multigrid method in connection with a Krylov solver substantially reduces the solution time of the electromagnetic problem in comparison to classical solution methods and can be effectively applied in parallel. Applications for the heat treatment of a gearwheel and for an automotive crankshaft are presented. The surface heat treatment of complex geometrical parts requires the introduction of a relative movement of workpiece and inductor to ensure a homogeneous surface treatment. A novel method is proposed, which is based on a discrete level set representation of the inductor motion that can be used to generate conforming finite element meshes in a Lagrangian setting.

Modélisation Numérique du

Keywords : Induction heating, surface hardening, Maxwell's equations, edge finite elements, level set method, algebraic multigrid

  Figure I.1 shows a typical CAD model of an automotive crankshaft. Two heat treated regions are indicated using arrows. The production process can be divided into pre-heating, focused induction heating, quenching and subsequently an optional tempering step.

Figure I. 1 :

 1 Figure I.1: Crankshaft with hardened features

Figure I. 2 :

 2 Figure I.2: Induction heating of a crankshaft 1 (Image provided courtesy of EFD)

Figure I. 5 :

 5 Figure I.5: Equilibrium diagram for non-alloyed low carbon steel C % Mn % Si % Ni % Cr % P % 0.41 0.71 0.22 0.26 1.04 0.022

Figure

  Figure I.6: Continuous Cooling Transformation diagram for 41Cr4

Figure I. 8 :

 8 Figure I.8: Current density per unit area for aluminium at 20 • C and f = 10 kHz

Figure

  Figure I.13: Hysteresis effect

Figure I. 15

 15 Figure I.15 visualises a first and second order phase transition. Figure I.15a shows the change of specific heat capacity during the change from solid to liquid phase of a typical one component system. The phase change happens under the condition of latent heat, which means that the phase transformation happens as isothermal process. It can be seen that C is not defined for point T pc . Standard engineering materials, like steel are multi-component systems. In these cases, the first order phase change does not happen as an isothermal process, but in a defined temperature range. Figure I.15b shows a second order phase change at the Curie

Figure I. 14 :

 14 Figure I.14: M and C vs. temperature T

  flux theorem relates the area integral of the electric flux density E with the volume integral of the electric charge density γ, enclosed by the former surface as

Figure II. 1 :

 1 Figure II.1: Two dimensional elliptic boundary value problem

Figure II. 2 :

 2 Figure II.2: Characteristics for parabolic PDEs

  Figure II.5: Isoparametric concept with reference triangle (a) and spatial triangle (b)

  . )dT e .

  introduces a base transformation ∂y ∂t = S ∂ ŷ ∂t and y = Sŷ, using matrix S, so that equation (II.122

Figure II. 6 :

 6 Figure II.6: Stability regions for the explicit Euler scheme (a) the implicit Euler scheme (b) and the trapezoidal rule (c)

  Figure II.7: Difference of sizes for the electromagnetic and heat diffusion problem

  .1. The streamlines of the magnetic field B for the first time step are indicated in figure II.8b.

Figure II. 8 :

 8 Figure II.8: Model of the ring inductor test case (a) and streamlines of the magnetic field (b)

  Figure II.9: Generation for a cutting surface on an example domain; initial finite element domain and bounded cutting plane (a), element cutting sets (b), final cutting surface (c)

Figure

  Figure II.10: Laplacian on a disc with hole; one cutting plane (a), gradient (b), detailed view of cutting planes and gradient (c), two cutting planes (d), gradient (e) and difference between gradients (f)

Figure

  Figure II.12: Example of the two cutting surface approach for computing conforming current sources; cutting surfaces (a), voltage potential (b) and streamlines of the source current (c)

Figure

  Figure II.13: Tree of the global domain (a), co-tree of the global domain (b) and top-view of the tree of the workpiece domain (c) for the ring inductor test case

  Figure II.14: Impact of full gauging on the convergence for the ring inductor test case using a two cutting surface based source current for an AMS, CLJP (a) and Jacobi (b) preconditioned GMRES solver

  Figure II.15: Computation of the voltage potential problem for the ring inductor using one virtual cutting plane (a), two virtual cutting planes (c) and resulting current distribution (b) and (d)

  Figure II.19 gives an overview of the post-processed magnetic flux, resulting from the first time step of the electromagnetic computa-

  Figure II.16: Convergence behaviour for different methods of including the source current

Figure

  Figure II.17: A scaled by the maximum value of case (f) for one cutting plane (ac), two cutting planes (d-f), with projection of source current on edge potential (b) and (d) and additional tree-gauging (c) and (f) using a Jacobi preconditioner

Figure

  Figure II.19: B scaled by the maximum value of case (f) for one cutting plane (ac), two cutting planes (d-f), with projection of source current on edge potential (b) and (d) and additional tree-gauging (c) and (f) using an AMS, CLJP preconditioner

Figure

  Figure III.1: Convergence behaviour for the solution of the magnetic vector potential formulation for the induction heating test case of a crankshaft with 3.5 million degrees of freedom; solver iterations for different preconditioned Krylov subspace methods (a); impact of different preconditioning methods for a GMRES solver (b)

Figure

  Figure III.2: multgrid discretisation on three levels (a) and sequence of operator applications for a standard V-cycle (b)

Figure

  Figure III.5: Setup time (a-e) and solver time (f-j) with respect to coarsening type for a partitioning on 4 processors

  Figure III.6: Setup time (a-e) and solver time (f-j) with respect to coarsening type for a partitioning on 16 processors

  Figure III.7: Convergence behaviour for different operator complexities for a partitioning on 4 processors

  Figure III.8: Convergence behaviour for different operator complexities for a partitioning on 16 processors

Figure

  Figure III.9: Impact of the operator ordering on the convergence behaviour

Figure

  Figure III.10: Distribution of discontinuous magnetic permeability µ r = f (B) of the crankshaft in the support region of the bearing, cut in the x-yplane

Figure

  Figure III.11: Distribution of discontinuous magnetic permeability µ r = f (B) of the crankshaft, cut in the y-z-plane

Figure

  Figure III.12: Solver iterations as a function of the relative magnetic permeability inside the crankshaft domain

Figure IV. 1 :

 1 Figure IV.1: Induction heating of a crankshaft (Image provided courtesy of EFD)

Figure IV. 2 :

 2 Figure IV.2: Model domain (a) and computational domain (b)

  t n+1 BV . Possibilities are bounding spheres, as well as simple bounding boxes. Figure IV.3b indicates such a bounding volume in green.

Figure

  Figure IV.3: Inital embedding domain (a) and spatially restricted search domain (b)

Figure IV. 4 :

 4 Figure IV.4: Restricted searchable finite element set (a) and identified removable element subset (b)

Figure IV. 5 :

 5 Figure IV.5: Recreated finite element mesh for the subsequent time step

  Figure IV.6a visualises the implicit level set function of a unit circle at coordinates x c = [1, 1], figure IV.6b gives that of a square with minimum point x min = [1, 1] and maximum point x max = [3, 3], whereas figure IV.6c shows a union of the former two implicit level set distance functions, which is defined as the minimal value.

Figure IV. 6 :

 6 Figure IV.6: Level set function of a circle (a), a square (b) and the union of these domains (c)

Figure IV. 7 :

 7 Figure IV.7: Sideview of a simple pyramid build from three tetrahedral elements.The third tetrahedron it is out of view.

  This section demonstrates the mesh adaptation algorithm for a heat treated automotive crankshaft. The location of the work piece and the inductor for initial time step t n = 0 s is shown in figure IV.8. The crankshaft is spatially fixed and depicted in blue. The inductor, which is shown in green, rotates around the local z-axis of the excentric part with an angular velocity of ω = 5 36 π rad s . The mesh is adapted for a time step of t n+1 = 0.1 s, which is equivalent to an inductor rotation of 2.5 • .

Figure

  Figure IV.8: Initial location of inductor and crankshaft for the mesh adaptation test case

Figure IV. 9 :

 9 Figure IV.9: Angle weighted pseudo normals for the inductor domain

Figure

  Figure IV.10: Closest point of a tetrahedral element with respect to the inductor domain

  Figure IV.11: Discrete level set function values for an inductor rotation of 2.5 • around the z-axis

Figure

  Figure IV.12: Cut of the recreated global finite element mesh after mesh adapation with new enclosing domain T t n+1 R

Figure

  Figure IV.13: Adapted finite element mesh T t n+1 corresponding to an inductor rotation of 2.5 • with temperature distribution of the crankshaft at time step t n+1 = 0.1 s

Figure V. 1 :

 1 Figure V.1: Model of the cylinder spin hardening test case

Figure V. 2 :

 2 Figure V.2: Temperature distribution (a), relative magnetic permeability (b) and heating power density scaled by the maximum value (c) for the cylindrical part at three different time steps

  Figure V.3: Streamlines of the magnetic field for the cylindrical part for three different time steps

  Figure V.4 shows some features of the geometry around which a circular ring shaped inductor is placed. The non-indicated corners have been rounded with a radius of 0.25 mm. The depth of this workpiece is 20 mm.

Figure V. 4 :

 4 Figure V.4: Geometrical details of the gearwheel test case

  (a) Location of the first cutting line (b) Location of the second cutting line

Figure V. 5 :

 5 Figure V.5: Finite element discretisation of the gearwheel geometry with indication of cutting lines for the heating energy density comparison

Figure V. 7 :

 7 Figure V.7: Temperature evolution for f = 1 kHz

Figure V. 8 :

 8 Figure V.8: Temperature evolution for f = 100 kHz

  Figure V.9: Convergence behaviour of the electromagnetic problem for the gearwheel spin hardening test case

Figure V. 10 :

 10 Figure V.10: Parallel convergence behaviour for the first electromagnetic time step of the gearwheel test case with f = 100 kHz for the CLJP (a)-(c) and HMIS (d)-(f) coarsening techniques

Figure V. 11 :

 11 Figure V.11: Crankshaft and inductor

Figure V. 12 :

 12 Figure V.12: Voltage potential of the inductor and placement of the flux concentrators

Figure V. 13 :

 13 Figure V.13: Cutting surfaces for the voltage potential computation of the closed inductor domain

Figure V. 14 :

 14 Figure V.14: Influence of flux concentrators on the distribution of the magnetic field for µ r = 1 for (a) and (c) and µ r = 1000 for (b) and (d)

Figure V. 15 :

 15 Figure V.15: Influence of flux concentrators on the distribution of the temperature for µ r = 1 (a), for µ r = 100 (b) and for µ r = 1000 (c)

  Figures V.20a and V.20b show the same region just after the Curie temperature has been reached in parts of the domain for time step ③, visible by the dark-blue shading in figure V.20b. The influence of the flux concentrators is clearly visible in the temperature evolution in figure V.20a, due to the large concentration of the high temperature directly below the flux concentrators.

Figure

  Figure V.17: Temperature (a), relative magnetic permeability (b) and electrical conductivity (c), shortly after the Curie temperature has been reached in parts of the domain

Figure V. 22 :

 22 Figure V.22: Comparison of induced heating power density for different distributions of relative magnetic permeability

  Figure V.23: Influence of the relative magnetic permeability of the flux concentrator on the convergence behaviour of the electromagnetic problem of the automotive crankshaft test case

  Figure V.24: Parallel convergence behaviour for the first electromagnetic time step of the crankshaft test case with µ r = 100 for the CLJP (a)-(c) and HMIS (d)-(f) coarsening techniques

  Le chauffage par induction électromagnétique est un procédé efficace permettant de chauffer directement une zone d'épaisseur contrôlée sous la surface de pièces métalliques en vue de les tremper. Cette thèse présente un modèle mathématique couplé électromagnétique/thermique et des approches numériques pour modéliser le procédé. Le modèle électromagnétique est basé sur une formulation en potentiel vecteur magnétique. Les courants de source sont imposés à l'aide d'une formulation en potentiel scalaire électrique permettant de modéliser des inducteurs de forme géométrique arbitraire. Le problème du transfert de chaleur est modélisé à l'aide de l'équation classique de diffusion de la chaleur. Le modèle électromagnétique est entièrement transitoire, afin de permettre l'introduction des effets non linéaires. La discrétisation spatiale est basée sur une approche éléments d'arêtes en utilisant un domaine global air/pièce/inducteur. Le système linéaire d'équations issu de la formulation implicite est creux et défini semi-positif ; il possède un noyau de taille importante. Il est démontré qu'un préconditionneur basé sur une méthode multigrille algébrique construit conjointement avec un solveur du type Krylov réduit substantiellement le temps de calcul du problème électromagnétique par rapport aux méthodes classiques de solution et peut être très efficace pour le calcul parallèle. Des exemples d'application pour le traitement thermique d'un pignon et pour un vilebrequin automobile sont présentés. Le traitement thermique des surfaces des pièces aux géométries complexes nécessite l'introduction d'un mouvement relatif de la pièce et de l'inducteur pour assurer un traitement homogène de la surface. Une nouvelle méthode est proposée, basée sur une représentation discrète d'une fonction level set du mouvement de l'inducteur qui peut être utilisée pour générer des maillages éléments finis conformes dans le cadre d'une configuration lagrangienne. Mots clés : Chauffage par induction, durcissement de surface, équations de Maxwell, éléments finis d'arêtes, méthode level set, multigrille algébrique Numerical Modelling of Induction Heating for Complex Geometrical Parts

  Le chapitre suivant introduit le modèle physique, décrivant les phénomènes électromagnétiques, fondé sur les équations de Maxwell, ainsi que le modèle classique de diffusion de chaleur. Ensuite, une version hyperbolique et une version parabolique sont dérivées pour les équations électromagnétiques. Conformément aux principes de la classification des équations aux dérivées partielles, en lien avec la méthode des caractéristiques, nous montrons pourquoi la version parabolique des équations électromagnétiques est plus adéquate pour décrire le problème du traitement par

induction. Ensuite, nous présentons brièvement la méthode de Galerkin en relation avec les éléments finis vectoriels du type Nédélec, ainsi que les formulations faibles pour chaque problème. L'aperçu sur des méthodes de discrétisation est finalisé par une introduction aux méthodes stables pour la discrétisation en temps. Finalement, un algorithme est présenté pour le calcul des courants conformes sur des géométries arbitraires et complexes, ce qui a été identifié comme condition indispensable pour que le problème électromagnétique soit bien posé.

  on ∂Ω D .

	(II.58)
	(II.59)
	(II.60)
	II.2.2.2 Weak form of the heat transfer equation
	The heat transfer equation, defined in its strong local form (II.34), can be refor-
	mulated into a weak form, like the elliptic voltage potential equation (II.30). The
	heat transfer problem is defined on the conducting domain Ω C . The second part of
	the left-hand side of equation (II.34) can be reformulated using Green's theorem,
	after multiplication by a test function v and integration by parts, such that

  Here, the conductor boundary is again divided into two separate subsets. Subset Ω CD defines the boundary with prescribed temperature T 0 , whereas Ω CN defines a boundary that permits heat exchange using a numerical flux, as specified in equations (II.35) and (II.36). The boundary domains are again distinct, such that∂Ω C = ∂Ω CN ∪ ∂Ω CD ,

		(II.62)
	with
	∂Ω CN ∩ ∂Ω CD = ∅.	(II.63)

  142)It can be noted that equation (II.142) is a projection on the curl of the complete test function space. Therefore, it is defined on the global domain Ω. Since no material quantities are present in equation (II.142), every finite element is qualitatively equal. It follows that standard preconditioned iterative solvers can be utilised to solve this equation. In this work, convergence could be rapidly achieved by using a simple Jacobian preconditioned conjugated gradient (CG) technique.

	Afterwards, the linear form b(v) of equation (II.80) can be replaced by the pro-
	jection
	b

Table II

 II

					∥ • 100%
		Jacobi AMS,CLJP Jacobi AMS,CLJP
	one cutting plane	407.98	407.99	698.89	698.88
	one cutting plane, Hrot	5.17	5.06	358.40	358.42
	one cutting plane, Hrot, tree	4.33	5.06	358.20	358.22
	two cutting planes	1.39	1.41	0.21	0.21
	two cutting planes, Hrot	1.08	6.10	0.005	0.002
	two cutting planes, Hrot, tree 0.0079	0	0.005	0

.2: Relative difference of the maximal values for the ring inductor test case

  sgn(a jk )

			a kj ,	else	(III.46)
	and
	âkm =		0, a km ,	if sgn(a km ) = sgn(a jm ) else.	(III.47)

  La fonction d'identification est basée sur une fonction discrète level set qui est créée en utilisant la forme primitive de l'inducteur discrétisé. L'algorithme est local, efficace et pourrait être parallélisé. Au début, le cadre du problème est défini. Ensuite, un aperçu est donné sur des méthodes existantes dans le contexte du traitement par induction et des méthodes venant d'autres applications. Les approches classiques, comme la méthode ALE ainsi que des algorithmes discrets sont présentées. L'avantage d'une approche avec la méthode locale de level set est montré. La description de l'algorithme est ensuite donnée. Un domaine de dimension deux est choisi pour illustrer la technique d'identification level set.

Ensuite, des résultats sont présentés pour un cas d'exemple : la rotation d'un inducteur autour d'un vilebrequin. Il est démontré que la qualité du maillage est préservée dans les étapes suivantes.

  Ce chapitre détaille quelques applications industrielles. Le but est de démontrer la fiabilité du modèle numérique qui a été présentée dans cet ouvrage. Les exemples sont rangés par ordre de difficulté croissante. Le premier exemple concerne un cylindre chauffé par un inducteur dans une configuration en boucle. Ce modèle est utilisé pour montrer l'application générale avec l'évolution du matériau en température. Les résultats sont comparés avec le modèle de pénétration exponentielle, qui a été introduit dans un chapitre précédent. Ensuite, le modèle numérique est appliqué pour analyser le traitement par induction d'un pignon. Il s'agit d'une géométrie complexe qui nécessite une discrétisation avec un maillage très raffiné. Le modèle est traité en supposant que les propriétés matériau restent constantes pendant la simulation. Néanmoins, il est montré pourquoi il serait nécessaire d'utiliser une simulation numérique pour cette géométrie complexe par rapport à une approximation simple en utilisant le modèle avec une estimation exponentielle. En effet, le courant d'annulation qui est créé dans les dents du pignon aboutit à une distribution complexe de la densité de chaleur générée par les courants induit.

	Le dernier cas présenté fait l'objet du projet OPTIPRO-INDUX et concerne le
	problème complètement non-linéaire du traitement par induction d'un vilebrequin.
	Nous analysons d'abord les d'application du courant sur un domaine complexe,
	i.e. l'inducteur utilisé. Nous analysons ensuite l'évolution de la température, avec
	prise en compte de l'évolution des propriétés matériau et une illustration des effets
	qui apparaissent. Finalement, l'efficacité numérique du solveur est illustré par sa
	convergence rapide et stable.

  The parallel convergence behaviour for this test case is indicated in figures V.10a, in connection with a CLJP coarsening technique, as well as for an HMIS convergence technique, in figure V.10d. It is visible that the HMIS convergence technique shows convergence problems in highly parallelised test cases. It can be noted that the convergence starts to severely degrade for parallel coarsening with 8 or more partitions.

	elements	nodes	edges	f	elec. time steps elec. solver time
	3.340.884 559.382 4.000.553 1 kHz	48	238 min
	3.340.884 559.382 4.000.553 100 kHz	48	412 min
	Table V.1: Discretisation, time steps and total cpu time for the electromagnetic
	computation of the gearwheel test case on one processor using a CLJP
	coarsening technique

  .2 for a sequential computation using a HMIS coarsening technique.Table V.2: Discretisation, time steps and total solver time for the fullly coupled electromagnetic/temperature evolution computation for the crankshaft test case with µ r = 100 on one processor using a HMIS coarsening technique

	elements	nodes	edges	electromagnetic time steps total solver time
	2.612.618 435.714 3.049.254	256	40.48h