Miren Cyprien

Luis

José Carlos

Victor, Geraldine Danae Antonio

Luca

José Alves

Edgar Pamela Simon

Keywords: Electromagnetic forming, high-speed forming, finite elements, image analysis

I would like to express my gratitude to my advisor professor François BAY. He gave me a lot of support in the theoretical, practical and also the human aspects of caring out a Ph.D research. He gave me many ideas, allowed me to knock on his door and discuss whenever I had a doubt. He helped me with many of the administrative aspects of being a foreign student in France. I also appreciated all the discussions we had on topics like politics, history, society and so on.

I could not leave aside Marie-Françoise GUENEGAN and Françoise TRUCAS

Dedicatory

A mi madre, Alice Xiomara Zapata Guevara Aun antes de haber empezado este sendero siempre soñaba con ilusión el día que escribiría las dedicatorias de mi tesis porque quería que pudieras leer estas lineas y así regalarte otro momento de orgullo y alegría.

Te fuiste justo antes de que pudiera dedicarte este detalle, este momento de felicidad. Aun así siempre llevo conmigo tu mirada de orgullo, tu mirada alentadora. Quiero decirte que este libro es tan mio como tuyo. Gracias a ti he recorrido el mundo, tu me enseñaste el valor de la independencia, a tenerme sobre mis propios pies.

Cada día de tu vida los dedicaste a mi hermana y a mi, hasta el ultimo momento continuaste tus luchas. Hoy ya encontraste reposo.

En tu honor, te quiero

To my mother, Alice Xiomara Zapata Guevara Before I began this path, I always dreamt of the day I would dedicate this work to you. I wanted you to read these lines and give you another moment of pride and happiness.

You left just before I could give you this little detail. Yet, I carry with me your gaze full of pride, full of encouragement. This book is as mine as it is yours. Thanks to you I have travel the world. You tough me the value of independence, to stand on my own feet.

Every day of your life you dedicated to my sister and me. Until the last moment you kept on fighting. Today you finally rest in peace.

iii On your honour, with love A mi padre, José Rodolfo Alves Contreras Dicen que el héroe de cada niño es su padre, hoy te confirmo que siempre has sido el mio. Desde pequeño admiraba tu inteligencia y sabiduría, para mi no había nadie persona mas inteligente en el mundo que mi padre. Tu me enseñaste a ser curioso, a amar el conocimiento pero sobre todo, a ser persistente y obstinado en la búsqueda de las respuestas.

Muchas personas a veces me preguntan como me pueden gustar tanto las ecuaciones y los números? Mi respuesta siempre es la misma: mi padre jugaba las matemáticas conmigo, él me enseño a verlas como algo divertido.

Este libro es para ti

Te quiero

To my father, José Rodolfo Alves Contreras Some people say that every boy's hero is his father. Today I want to confirm that you have always been mine. Since I was a child, I admired your intelligence, your wisdom, for me there was no one as intelligent as my father. You tough me to be curious, to love knowledge, to be persistent in the research for answers.

Many people ask me how can I like so much equations and numbers ? My answer is always the same: my father tough me math as a game, teaching me to see it as something fun.

This book is for you, With love iv

Résumé: Introduction

Nous introduisons ici les principaux enjeux associés aux procédés de formage par impulsion magnétique. Nous essayons de répondre à deux questions fondamentales:

• pourquoi? Ou l'intérêt industriel d'un tel procédé

• comment cela fonctionne.

Nous présentons en parallèle un historique du procédé. Nous concluons ce chapitre par les principaux objectifs de ce travail de thèse.

Le formage électromagnétique est une technologie utilisée de façon corissante dans les dernières décennies -grâce notamment à la formabilité accrue qu'il offre pour les matériaux à haute résistance et faible masse spécifique tels que les alliages d'aluminium et de magnésium.

Ce procédé rentre dans la catégorie des procédés de formage à haute vitesse de déformation. La déformation de la matière se produit dans une plage entre 103 et 104 s-1. Plusieurs études montrent que l'on obtient dans cette gamme de vitesses des limites de formabilité bien plus élevées que celles rencontrées dans le cadre d'un procédé quasi-statique.

Un autre avantage important réside dans l'élimination presque totale du retour élastique. Ce phénomène est en effet un défi récurrent dans la mise en forme des pièces dès lors que nous utilisons des procédés basés sur la déformation plastique.

Le procédé a vu le jour aux alentours des années 1920 avec les études de Kapitza Kapitza 1924. L'industrialisation du procédé a du attendre les années 50-60, lorsque, suite à un accident dans une installation nucléaire aux États-Unis on a trouvé le moyen de générer les niveaux d'énergie requis pour déformer la matière à haute vitesse.

Depuis, plusieurs études ont été accomplies pour essayer de comprendre comment l'interaction entre les différents aspects physiques conduit aux résultats obtenus dans la pratique.

Le présent travail consiste principalement dans le développement d'un modèle numérique basé sur la méthode des éléments finis qui permet de prendre en compte les différents aspects multiphysiques dans un cadre complètement tridimensionnel. Nous allons développer un modèle qui soit à la fois précis pour faire une étude scientifique approfondie et en même temps avec une interface "simple" pour permettre l'utilisation efficace de l'outil numérique dans le cadre industriel de conception d'un procédé.

Chapter 1. Introduction

Interest in lightweight structural applications has driven research on alternative forming processes that allows costs reduction while maintaining or improving structural resistance. The first stage for introducing lightweight components into the fabrication line is shape optimization because structural resistance can be highly upgraded by optimally adapting the geometry of the component to the application. This optimal shape leads in many cases to complex geometries that cannot be easily obtained by traditional forming processes.

The second stage consists in selecting materials which satisfy both the weight requirements and the mechanical resistance requirements. In structural applications for the aeronautic and automotive industries the most popular options are aluminium and magnesium alloys. Nevertheless, most aluminium alloys show poor formability qualities when compared to steel using classical (low speed) processes. [Daehn 2002] enumerates some of the main difficulties associated with the forming of this material:

1. Forming limits of aluminium are significantly lower than those for steel. Aluminium is particularly prone to tearing at bends. This limits the shapes that can be fabricated and slows die design, die tryout and application.

2. As the elastic modulus of aluminium is lower than that of steel, springback is more severe and it is difficult to keep dimensional tolerances.

3. Conventional die try-out with mating male and female dies is always slow and expensive. It would be desirable if a least prototypes could be produced with one-sided die sets. This would speed time-to-market.

Several experimental works show that it is possible to access increased formability ranges and almost eliminate the springback by means of high speed forming technologies. These qualities have already been used in several joining and forming applications as shown in fig. 1.1. Besides, the main high speed technologies used up to now require one-sided die due to the intrinsic way of delivering the load.

• Explosive forming: A chemical or combustible charge is placed near the work piece inside a confined chamber. When the load is ignited the pressure shock wave transmitted through the air will cause the high speed deformation of the part. • Electro-Hydraulic forming (EHF): An electrical arc is formed between two conductors placed near each other inside a container with water.

The pulsed nature of the arc will induce a fast phase transformation from liquid to high temperature vapour. The high energy shock wave produced by this phenomenon is responsible for the forming process.

• Electromagnetic forming (EMF): Induced Lorentz forces on the work part due to a fast electric discharge into a coil placed near the piece lead to a high acceleration and the resulting deformation of the material. The basic principle behind all these processes is, as depicted in fig. 1.2, that given a certain amount of energy (stored as chemical or electrical energy), the delivered power is increased by reducing the discharge time. Consequently, Chapter 1. Introduction the loads generated on the material are very high without the need of larger1 amounts of energy.

Important requirements for industrializing a given technology are its controllability and predictability. Safety and pollution are also major factors that will help in the choice of the most appropriate and "friendly" technology. This has driven attention towards EHF and EMF because they offer access to the same range of high speed loads as explosive forming, but with increased controllability. Besides, they are much more environmental friendly (no chemical reactive and no need for lubrication).

In this work we have chosen to focus our attention on the last of these three technologies, electromagnetic forming. From a scientific point of view we are interested in studying the interaction between the different physics. Full understanding of the process is quite challenging because of the multiple physical aspects involved: namely, the electromagnetic wave propagation, fast induction heating and self-heating through deformation work, the inertial effects on the formability, etc. And from a technological point of view, we want to develop a numerical tool that provides accurate representation of the complete process.

The magnetic pulse forming process

The magnetic pulse forming (MPF) technology or electromagnetic forming (EMF) consists in deforming metallic components through the application of an intense electromagnetic pulse. This process enters in the category of high speed forming processes because of the range of strain rates that are usually attained, ranging from 10 3 s -1 to 10 4 s -1 .

It has been one of the most commercially used high speed processes for the last 40 years because it offers high controllability and safety, making it apt for industrialization and automation.

Motivation and Challenges

Motivation and Challenges

Some of the most attractive assets from this technology have been identified by [START_REF] Psyk | Electromagnetic forming -A review[END_REF]. Here, we summarize some of these points:

• It is a one-sided die process. The coil occupies the roll of the "male die". Furthermore, a single coil can be used for several final geometries.

• Since the forces are applied without any mechanical contact it is possible to perform surface treatments or rectifications before the forming process is applied. In metal-to-metal welding applications it is not necessary to "clean" the contact surfaces. During impact, a jet stream is generated in the collision zone that cleans the surfaces (fig. 1.2), removing the oxides and thus creating a high quality bonding BMAX 2014. • It is possible to join dissimilar materials that could not normally be joined by classical welding. For metals this is achieved by the atomic bonding due to the high speed impact. Joining without welding is possible by the almost perfect incrustation of the work piece into the receptor. At the high strain rates of the process the material undergoes a transition from elastic-plastic behaviour to a more hydrodynamic behaviour.

• One of the most interesting yet poorly understood outcomes is the reduction (almost full elimination) of the elastic spring back -enabling die design with the exact final shape.

• It can be tagged as a highly environmental friendly process because the lack of moving parts and the lack of contact with the work piece eliminates the need for lubrication.

• It is also considered as a low cost process. Assuming that a single "shot" consumes 10kW from the grid. The charging time of the capacitor being Chapter 1. Introduction around 5s, and using the reference price of 0.144 euros /kW h (July 2015, EDF, France). We obtain an approximate price of 0.002 euros per shot.

As any other technology, EMF has its limitations or disadvantages:

• It is limited to use materials with high electric conductivity. One workaround is the use of a punch or driver. A conductive material is placed near the inductor to which the kinetic energy is given. The punch is accelerated and forms the low-conductive part by impact. Another way is to increase the discharge frequency of the capacitor in order to diminish the skin depth penetration of the fields in the work part, thus avoiding the fields to traverse and create forces on the opposite side that will work against the acceleration of the material.

• The coil has to be carefully designed because it will withstand large thermal and mechanical loads. Given a certain frequency of production, fatigue can be a serious damaging factor. This is one of the most limiting factors of the process.

• While it was shown before that the process is profitable from an economical point of view, it is also quite ineffective in terms of energy recovery. It is reported that the rate of effective plastic energy used to the stored energy in the capacitors is between [2 -20]%, normally closer to 2%. One of the causes being the resistive losses in the electric connections. Another cause could be the energy elastic losses involved at impact.

History

The first time someone effectively used such a process to deform a solid conductor beyond its yield point was Kapitza 1924. He was able to produce electromagnetic fields about 15 to 20 kiloGauss (compared to the earth's natural magnetic field of 0.2Gauss) but the technological constraints at that time required enormous power sources for the generation of the energy and also for the cooling of the inductors which was already a major blocking point. So the technology was abandoned for almost 40 years.

As explained in Zittel 2010 it was at the end of the 1950's that during a nuclear fusion research by General Atomic2 in San Diego, material failures were de-

History

tected between the conductor carrying the electric current. This accident led to the patent of [START_REF] Harvey | Metal forming device and method[END_REF]. This marked the start point for the industrialization of the process. In this work, they already stated the main classical forming configurations (mainly related to the geometry of the coil) that can be found in industry, see fig. 1.4. We can find tube forming by either compression or expansion, in which case a cylindrical coil is needed. We also have flat sheet forming or bulging. In this case a flat or "pancake-like" coil is used. The technique of indirect forming within the EMF process was published by Brower 1966 in which an elastomeric punch is used to transmit the energy. Thus allowing to form non-conductive materials. Some of the first work on joining of dissimilar materials were done by [START_REF] Al-Hassani | Analysis of electromagnetic forming process[END_REF] to join a metallic cap with a glass bottle and, [START_REF] Rafailoff | Fügen von rohrförmigen Werkstücken aus Metall und Porzellan durch Magnetumformung[END_REF] where a metallic tube was joined with a porcelain component. A more detailed account on the evolution of the technology can be found in [START_REF] Psyk | Electromagnetic forming -A review[END_REF].

As the technology gained interest in the scientific and industrial domains, the need for accurate and representative tools for predicting the outcomes increased. According to Daehn 2002 between the end of the 60s and the 70s most of the "modelling" was done in closed-forms (0-D approximations) enabling to compute scalar quantities, basically estimations of the order of magnitudes. One of the works documented here is the 1977 Ukranian Electromagnetic Forming Handbook. In the same line, the work of [START_REF] Leroy | Formage électromagnétique[END_REF] served as an abaqus in the 1980's in French industry for several industrial applications.

At the 1-D level we find the equivalent circuit approaches which are also interesting, mainly for the computation of the electric system. In this case only simple geometries can be treated, the EM fields are treated in average terms and the mechanical deformation phenomena is not considered.

Chapter 1. Introduction

One of the first attempts to compare experimental results with 2-D finite elements/finite difference computations was done by [START_REF] Takatsu | High-speed forming of metal sheets by electromagnetic force[END_REF]. They studied a free bulging case of an aluminium sheet in order to find a constitutive law for the plastic behaviour of the material. This case has been largely referenced thereafter.

A few years later appeared one of the first documented 2-D codes dedicated to EMF. It was the CALE code (C Language Arbitrary Lagrangian-Eulerian), mainly dedicated for research activities. In the work of Fenton 1996 they successfully applied this method for the computation of axi-symmetric models in order to study plasticity in the high strain rates range.

In recent years, there have been several works in the development of 3-D tools for the simulations of the process. The two main toolbox found in literature are ANSYS (see Siddiqui 2009;[START_REF] Bartels | Simulation Models of the Electromagnetic Forming Process[END_REF][START_REF] Demir | Simulation of Wrinkle Formation in Free Electromagnetic Tube Compression[END_REF][START_REF] Cui | 3D modeling and deformation analysis for electromagnetic sheet forming process[END_REF] for instance) and LS-DYNA (L' Eplattenier et al. 2009;[START_REF] L'eplattenier | Update on the Electromagnetism Module in LS-DYNA". 12th LS-DYNA Users Conference[END_REF]. Nevertheless, industrial use of the 3D tools have had a slow introduction in industry due to the complexity of setting up the simulations and the computational time required in order to accomplish a coupled simulation.

Understanding how it works

In order to design a target oriented toolbox for the process, we need to understand how it all works. The process can be decomposed in the following consecutive stages: 1. The stored energy in the capacitor bank is rapidly discharged as a pulse signal into the the system.

1.5. On the increased formability properties 2. The time-varying current in the coil generates intense electromagnetic fields near the work piece.

3. Induced eddy currents arise on the work piece flowing in opposite direction to the current in the coil.

4. As a consequence, repulsive Lorentz forces appear on both parts.

5. The forming takes place by the transformation of the electromagnetic force into kinetic energy, and subsequently into plastic deformation at high strain rates. In joining applications, the high speed impact is the responsible of bonding the surfaces at the atomic level.

On the increased formability properties

Several works on high speed forming of certain metals alloys (aluminium and copper mainly) have shown that for the given loading conditions increased formability properties are found as seen in fig. 1.6. In [START_REF] Balanethiram | Hyperplasticity: increased forming limits at high workpiece velocity[END_REF] two possible causes are enumerated: (1) A change in the material constitutive law that makes the material more formable and (2) The effect of inertia attenuating the growth of cracks. It has been noted that the most relevant factor is the inertia that stabilizes the development of necking and increases the global ductility. In the works of Chapter 1. Introduction [START_REF] Jacques | Effects of microscale inertia on dynamic ductile crack growth[END_REF][START_REF] Molinari | Dynamic failure of ductile materials[END_REF] a model of micro-inertia effects was developed in order to study the influence of inertia on the propagation of cracks. It was concluded that the micro-inertia helps reducing voids growth (hence the coalescence) in the regions where strain is more concentrated, thus allowing a more progressive accumulation of damage, which is distributed in several regions rather than concentrating in the most fragile zone as would be expected under low speed loading conditions.

Objectives

Following the observations made above it was detected that there is a need for developing a simulation toolbox that is capable of handle the multiphysical complexities linked to the process, and at the same time be highly predictive.

In order to accomplish such task we subdivided the study in the following parts:

1. Include the resolution of the full Maxwell equations into the electromagnetic module MATELEC. The solver was originally conceived for induction heating applications, meaning applications in a low frequency regime for which several hypothesis are made that do not match the needs for EMF processes.

2. Include the EMF machine into the simulation set-up using an electric circuit description that is coupled with the finite elements resolution. The methodology behind these two points are described in chapter 3.

3. Optimize the simulation CPU time by means of parallel resolution techniques. For this purpose we have selected to focus the simulation strategy of the electromagnetic problem on full finite element discretization in contrast to boundary elements/finite elements couplings.

4. Couple the thermo-mechanical (TM) solver with the EM solver in order to guarantee an efficient communication and accurate reproduction of the phenomena within a reasonable CPU time. The numerical resolution of the physics are formulated using implicit methods in order to guarantee the numerical accuracy3 .

1.6. Objectives 5. One of the difficulties in the simulation of the electromagnetic phenomena within the process is the management of movement/deformation. chapter 4 will provide the details of the adopted methodology.

6. Compare the simulation results with available data in literature and also with experiments performed in our facilities.

Résumé: Modèles physiques et mathématiques

Ce chapitre est consacré à la description mathématique des différentes problématiques physiques impliquées: la mécanique des solides, le transfert de chaleur et la propagation des ondes électromagnétiques dans le cadre macroscopique de la mécanique des milieux continus.

Nous démarrons avec l'introduction d'un schéma d'interaction entre les trois phénomènes physiques évoqués. Nous utilisons ensuite les concepts de la mécanique des milieux continus pour développer les équations différentielles fondamentales à partir des équations intégrales de conservation.

Une fois établies les équations différentielles de conservation, nous rappelons certains modèles physiques pour les équations constitutives des matériaux, qui permettent de lier les différentes quantités physiques et qui permettent donc de décrire le comportement de la matière.

A coupled multiphysics problem

The electromagnetic forming process can be decomposed in three main elements (fig. 2.1):

• The EMF machine, consisting in the command bank and the electric capacitors bank for energy storage. This element contains no moving parts since its purpose is the set-up of the process and energy delivery through electric connections. From this fact it is possible to reduce the modelling of it by only considering the electric parameters.

• The work bench: Its main component is the coil, responsible of transmitting the electromagnetic waves. In case of forming applications we will also have the matrix or female-die.

• The work piece, which is the target element, subject to the deformation by either direct or indirect forming. For the simulation of the whole process we shall consider the interaction of these different elements. In order to accomplish this task we will review the mathematical description behind the physics governing the process.

A coupled multiphysics problem

The core complexity related to this process lies in its multiphysical nature. All the different physical phenomena are highly interconnected. The electromagnetism drives the process by creating the "electromagnetic pressure" on the work piece. As the material deforms and changes its temperature due to the Joule heating effect and plastic self-heating, the material properties change including the electromagnetic ones. Also it is important to closely keep track of Before establishing the schemes for the interaction between the physics we are going to review the mathematical notation used to describe the different problems. The coupling strategy will be discussed in the following chapters together with the finite elements approach.

Continuum Mathematics

This section will present some of the basic notations and concepts related to the mathematics used in continuum mechanics for the further development of the equations.

Reference frames

The description of the continuum mechanics phenomena starts by defining the reference systems shown in fig. 2.3. The Eulerian or laboratory frame allows measuring the quantities with tools considered static with respect to the observer. The Lagrangian frame follows the movement and deformation of the given domain Ω in space so as to imagine that for an observer in such a frame, no motion has taken place.

Under this framework, the kinematics of a selected point P ∈ Ω with coordinates #» x = #» x ( #» x , t) are described by Displacement

Continuum Mathematics

#» u = #» u ( #» x , t) = #» x ( #» x , t) -#» x ( #» x , t 0 ) = #» x ( #» x , t) -#» x (2.1) Velocity #» v = d dt #» x ( #» x , t) = d dt [ #» u + #» x ] = d #» u dt = ∂ #» u ∂t + ∂ #» u ∂ #» x d #» x dt = ∂ #» u ∂t (2.2) Acceleration #» a = d dt #» v = ∂ #» v ∂t + ( #» v • ∇) #» v = ∂ ∂t + #» v • ∇ #» v (2.3)
The measure equivalences between the deformed and undeformed configurations are given by:

Length d #» x = dx i = ∂x i ∂x j dx j = F ij dx j = F d #» x (2.4) Volume dV = |dx êx • (dy êy × dz êz )| = dx dy dz (2.5) dv = |dxê x • (dyê y × dzê z )| = det F dV = J dV (2.6) Surface ndS = J F -T n dS (2.7) (Nanson's Formula)
Other additional concepts that are going to be necessary are 1

1 See demonstration in section A.1 for the derivative of the Jacobian

Chapter 2. Physical & Mathematical Framework Velocity Gradient Tensor d #» v = ∂ ∂ #» x ∂ #» x ∂t • ∂ #» x ∂ #» x = ∂ ∂t ∂ #» x ∂ #» x • ∂ #» x ∂ #» x (2.8) d #» v = ∂ ∂t F • F -1 • d #» x = Ḟ • F -1 d #» x = L • d #» x L = Ḟ • F -1 = ∂ #» v ∂ #» x
Derivative of the Jacobian

d dt J = ∂ ∂t J = J Tr F -1 • Ḟ = J ∇ • #» v (2.9) Theorem 1 (Theorem of conservation). Given a functional F existing over a domain Ω, if F is C 1 and if/∀v ⊂ Ω v F dΩ = 0 ⇒ F = 0 ∈ Ω Theorem 2 (Gauss Theorem). Given a function f : R 3 → R 3 , if f ∈ C 1 Ω ∇ • f dΩ = Γ f • n dS (2.10)
In case of a second order tensor we have

Ω ∇ • T dΩ = Γ T • n dS (2.11) Theorem 3 (Stokes' Theorem). C #» F • d #» l = Γ ∇ × #» F • n dS (2.12)

Reynold's Transport theorem -Volume Integral time derivative

As will be shown later, many of the fundamental equations of mechanics are obtained by establishing the conservation laws and then deriving the differential forms from integral representation.

The Electromagnetic Problem

Derivative of a volume scalar integral

d dt Ω f dΩ = lim ∆t→0 1 ∆t    Ω(τ +∆t) f ( #» x , τ + ∆t) dΩ - Ω(τ ) f ( #» x , τ ) dΩ    = lim ∆t→0 1 ∆t   Ω 0 f ( #» x , τ + ∆t)J ( #» x , τ + ∆t) dΩ - Ω 0 f ( #» x , τ )J ( #» x , τ ) dΩ  
Now, the integration limits are time independent, allowing to write

d dt Ω f dΩ = Ω 0 ∂ ∂t (f ( #» x , τ )J ( #» x , τ )) dΩ = Ω 0 ∂f ∂t J + f ∂J ∂t dΩ = Ω 0 ∂f ∂t J + f J ∇ #» v dΩ = Ω 0 ∂f ∂t + f ∇ #» v J dΩ = Ω df dt + f ∇ #» v dΩ = Ω ∂f ∂t + #» v • ∇f + f ∇ #» v dΩ = Ω ∂f ∂t + ∇( #» v f ) dΩ (2.13) d dt Ω f dΩ = Ω ∂f ∂t dΩ + Γ f #» v • n dS (2.14)
From (2.14) we find Reynolds' transport theorem.

The Electromagnetic Problem

Note: In this section we are going to perform a non-classical derivation of the equations of electromagnetism following what is known as the Heaviside-Hertz model. The objective being to point out that there are several questions regarding the development of Maxwell equations that are left without answer or are overlooked in most textbooks. More than a matter of being wrong or not, we rise the question on the mathematical steps for the derivation of the model. After realizing that Maxwell equations do represent the electromag-Chapter 2. Physical & Mathematical Framework netic phenomena, modern physics has taken the position of saying that "it does not matters how the model is obtained, what matters are the equations themselves". If we revisit the equations from the beginning, could we find some harmony in the procedure? A more "logical" derivation?... We will see that using the Heaviside-Hertz model we can follow the tools from continuum mechanics to derive the complete set of equations. It is also found that the use of this model does not violate the principle of special relativity!. Nevertheless, since this topic goes beyond the objectives of the present work we shall simply point out our findings, doubts and possible path for the discussion. The numerical development of the rest of the work will still use the Maxwell model given that at the work conditions, both models are equivalent.

The critics to the Maxwell model are summarized here:

1. Before Maxwell, it was accepted that the law for the magnetic field produced by steady currents was

∇ × #» B = #» J ε 0 C 2
Maxwell realized that there was a problem with such representation. It violates the principle of charge conservation [START_REF] Feynman | The Feynman lectures on physics[END_REF].

∇ • #» J = - ∂ρ e ∂t
So he proposed to add the term ∂ #» E ∂t to the r.h.s of the equation to get

c 2 ∇ × #» B = #» J ε 0 + ∂ #» E ∂t
With this addition Maxwell had fully represented the electromagnetic phenomena. However, Einstein 1905 in his paper "On the Electrodynamics of moving bodies" starts by saying It is known that Maxwell's electrodynamics -as usually understood at the present time-when applied to moving bodies, leads to asymmetries which do not appear to be inherent in the phenomena.(...) Making reference to the observed phenomena between a magnet and a conductor whilst moving one relative to the other. This observation in addition to the constancy of the velocity of light led to the development of the special relativity.

The Electromagnetic Problem

2. It has been criticized that the force on a charge particle can not be directly obtained from the solution of the field equations but has to be postulated as the "fifth" equations of electrodynamics.

#» F = q( #» E + #» v × #» B)
The development that follows will mainly tackle both critics presented before.

Electromagnetic Equations

Let's start from the integral basic laws Faraday's Induction Law

C #» E • d #» l = - d dt Γ #» B • n dS (2.15a) Maxwell-Ampere's Law C #» H • d #» l = - d dt Γ #» D • n dS (2.15b) Gauss's Law Γ #» D • n = Ω ρ e dΩ (2.15c)
Gauss's Law for the magnetic field Following Maxwell's model we would have done

d dt Γ #» B • n dS = ∂ ∂t Γ #» B • n dS = Γ ∂ #» B ∂t • n dS (2.16) Chapter 2. Physical & Mathematical Framework
Which can be done if the domain of integration is independent of time, namely, if it doesn't move. Let us instead follow the procedure for the time derivative of a volume integral. Assume #» A( #» x , t) be a given vector field, then

d dt Γ #» A • n dS = lim ∆t→0 1 ∆t    Γ(τ +∆t) #» A( #» x , τ + ∆t) • n dS - Γ(τ ) #» A( #» x , τ ) • n dS    = lim ∆t→0 1 ∆t   Γ 0 J F -T #» A( #» X, τ + ∆t) • n dS 0 - Γ 0 J F -T #» A( #» X, τ ) • n dS 0   = Γ 0 nT • ∂ ∂t J F -1 #» A( #» X, τ ) dS 0 = Γ 0 nT • J F -1 ∂ ∂t #» A + ∂J ∂t F -1 #» A + J ∂ ∂t F -1 #» A dS 0 = Γ 0 nT • J F -1 ∂ ∂t #» A + (J ∇ • #» v )F -1 #» A -J F -1 Ḟ F -1 #» A dS 0 = Γ 0 J F -T n • ∂ ∂t #» A + (∇ • #» v ) #» A -∇ #» v • #» A dS 0 (2.17) = Γ n • d dt #» A + (∇ • #» v ) #» A -∇ #» v • #» A dS = Γ n • d dt #» A + (∇ • #» v ) #» A -( #» A • ∇) #» v dS (see 2 ) = Γ n • ∂ ∂t #» A + ( #» v • ∇) #» A + (∇ • #» v ) #» A -( #» A • ∇) #» v dS = Γ n • D Dt #» A dS (2.18)
From (2.18) we retrieve an equivalent of the Upper Convected Oldroyd derivative for a vector field.

D Dt

#» A = ∂ ∂t #» A + ( #» v • ∇) #» A + (∇ • #» v ) #» A -( #» A • ∇) #» v (2.19)
In Christov 2006 was implied that by choosing this time derivative definition, it was possible to achieve the time invariance of the Hertz's equations. Here,

2 ∇ #» v • #» A = ∂ i (v j êj ) ⊗ êi • A k êk = A k ∂ i v j êj ⊗ êi • êk = A k ∂ i v j êj δ ik = A i ∂ i v j êj = ( #» A • ∇) #» v

The Electromagnetic Problem

we stress out that (1) More than a matter of choice, it is the natural outcome by construction (2) Time invariance is not a consequence of the development, from our perspective, it is just a misinterpretation based on the unconscious assumption of F = I.

By means of the vector calculus identity

∇ × ( #» A × #» v ) = ( #» v • ∇) #» A + #» A(∇ • #» v ) -( #» A • ∇) #» v -#» v (∇ • #» A) (2.20)
we can reduce (2.19) as

D Dt #» A = ∂ ∂t #» A -∇ × ( #» v × #» A) + #» v (∇ • #» A) (2.
∇ × #» E = - ∂ #» B ∂t + ∇ × ( #» v × #» B) (2.22a) ∇ × #» H = ∂ #» D ∂t -∇ × ( #» v × #» D) + ρ e #» v (2.22b) ∇ • #» D = ρ e (2.22c) ∇ • #» B = 0 (2.22d)
In (2.22) we find the electric current density term

#» J = ρ e #» v
Another fundamental equation is the equation of charge continuity in a moving frame of reference. It can be obtained simply by taking the divergence of (2.22b)

∇ • ∇ × #» H = ∇ • ∂ #» D ∂t -∇ × ( #» v × #» D) + #» J 0 = ∂ ∂t ρ e + ∇ • ( #» J ) (2.23) Chapter 2. Physical & Mathematical Framework

Electromagnetic Force -The Lorentz Force

As argued before, it has been criticized that the Lorentz force cannot be obtained as a solution of the field equations, but has to be postulated separately.

In [Phipps 2004] it is shown that the Lorentz force can be derived from the Hertz field equations without any additional postulate. An important remark is that in their derivation, they find an additional constant value, for which several arguments concluded that its value should be zero. The cause is that instead of starting from the integral equations, they simply took the Maxwell equations and replaced

∂ ∂t → d dt = ∂ ∂t + ( #» v • ∇).
Let's instead simply take the definition of the electric field, hence establish that for a point charge q (or a distribution of charges ρ e )

#» F L = q #» E or #» f L = ρ e #» E (2.24)
It is possible to demonstrate that this force is the same as the one given by the Maxwell representation using the potential formulation.

Potential Formulation An usual way to solve the differential equations is through the definition of the scalar and the vector potentials. This methodology can also be applied to (2.22), finding

∇ • #» B = 0 ⇒ #» B = ∇ × #» A (2.25) ∇ × #» E = - ∂ #» B ∂t + ∇ × ( #» v × #» B) ⇒ ∇ × #» E = - ∂ ∂t (∇ × #» A) + ∇ × ( #» v × ∇ × #» A) ⇒ ∇ × #» E + ∂ #» A ∂t -#» v × ∇ × #» A = 0 (2.26)
Since for any scalar function φ, ∇ × (-∇φ) = 0 holds, then

⇒ #» E + ∂ #» A ∂t -#» v × ∇ × #» A = -∇φ ⇒ #» E = -∇φ - ∂ #» A ∂t + #» v × ∇ × #» A (2.27)
Where, for the magnetic field, we find the traditional definition

#» B = ∇ × #» A but

The Electromagnetic Problem

for the electric field there is an additional term

#» v × ∇ × #» A.
We thus find that using (2.27)

#» F = q -∇φ - ∂ #» A ∂t + #» v × ∇ × #» A #» F = q #» E M axwell + #» v × #» B (2.28)
The "M axwell" subscript is meant to highlight that the fields are in fact different quantities in each theory. But the effect over matter remains unchanged (the force) if we use the potential description as reference.

Classical electromagnetism

The classical representation of the electromagnetic problem, as given by Maxwell can be found by using the integral representation of the Maxwell-Ampere's equation

C #» H • d #» l = - d dt Γ #» D • n dS + Γ #» J • n dS (2.29)
Then, the l.h.s. is modified using Stokes theorem and the r.h.s by saying

d dt Γ #» D• n dS = Γ ∂ ∂t #» D • n dS.
This results in a differential set similar to eq. (2.22) where the difference is that the terms related to the velocity are gone. This results in

∇ × #» E = - ∂ #» B ∂t (2.30a) ∇ × #» H = ∂ #» D ∂t + #» J (2.30b) ∇ • #» D = ρ e (2.30c) ∇ • #» B = 0 (2.30d) Chapter 2. Physical & Mathematical Framework

The Solid Mechanics and Heat Transfer Problems

At the macroscopic scale, the fundamental laws for continuum media is given by the following four conservation laws

Mass conservation:

The mass contained in a given domain Ω is given by m = Ω ρ( #» x , t) dΩ, where ρ is the density of matter at each point of the domain. Assuming no flow of matter on the boundaries of the domain (and neglecting the transformation of mass into energy)

Dm Dt = D Dt Ω ρ( #» x , t) dΩ = 0 (2.31)
In a Lagrangian description we can write

Ω ρ dΩ = Ω 0 ρ 0 dΩ ⇒ ρJ = ρ 0 , J = det F (2.32)
In the Eulerian frame we obtain the conservative form

∂ρ ∂t + ∇(ρ #» v ) = 0 (2.33)
Equation (2.33) is used to define the compressibility condition which is generally expressed as

∂ρ ∂t + #» v • ∇ρ + ρ∇ • #» v = 0 (2.34) dρ dt + ρ∇ • #» v = 0 ∇ • #» v = - 1 ρ dρ dt (2.35)

Linear momentum conservation:

If we now consider that a set of volumetric forces (ρ #» b ) and superficial tensions ( #» t ) act upon the domain Ω, the total forces are given by 2.4. The Solid Mechanics and Heat Transfer Problems

f (t) = Ω ρ #» b dΩ + Γ #» t dS (2.36)
The total linear momentum of the body is

#» P (t) = Ω ρ #» v dΩ (2.37)
From Newton's second law

D #» P (t) Dt = f (t) ⇒ D Dt Ω ρ #» v dΩ = Ω ρ #» b dΩ + Γ #» t dS (2.38)
The r.h.s of (2.38) is given in two forms, each one more convenient for different applications

D Dt Ω ρ #» v dΩ = Ω ρ D #» v Dt dΩ (2.39a) = Ω ∂ ∂t (ρ #» v ) + ∇(ρ #» v • #» v ) dΩ (2.39b)
Superficial tension vector #» t is measured in a plane defined by the normal n is given in term of the Cauchy stress tensor as

#» t = σ • n (2.40)
Then, by means of the Reynolds' transport theorem, the l.h.s is

Ω ρ #» b dΩ + Γ #» t dS = Ω ρ #» b dΩ + Γ σ • n dS = Ω ρ #» b dΩ + Ω ∇ • σ dΩ
The stress tensor is usually decomposed in a deviatoric and a hydrostatic com- From where

Ω ρ #» b dΩ + Γ #» t dS = Ω ρ #» b + ∇ • S -∇p dΩ (2.43)
We thus obtain

Lagrangian ρ D #» v Dt = ∇ • S -∇p + ρ #» b (2.44) Eulerian ∂ ∂t (ρ #» v ) + ∇(ρ #» v • #» v ) = ∇ • S -∇p + ρ #» b (2.45)

Angular momentum conservation:

The angular momentum ( #» h ) of the domain is given by

#» h (t) = Ω #» x × ρ #» v dΩ (2.46)
The total momentum from the external loads is defined as

# » M (t) = Ω #» x × ρ #» b dΩ + Γ #» x × #» t dS (2.47)

The Solid Mechanics and Heat Transfer Problems

Then, the conservation law says

D #» h (t) Dt = # » M (t) D Dt Ω #» x × ρ #» v dΩ = Ω #» x × ρ #» b dΩ + Γ #» x × #» t dS Ω ρ D Dt ( #» x × #» v ) dΩ = Ω #» x × ρ #» b + #» x × ∇ • σ + : σ dΩ ⇒ Ω #» x × ρ D #» v Dt -∇ • σ -ρ #» b Linear Momentum dΩ = 0 = Ω : σ dΩ (2.48)
Where is the Levi-Civita operator. From here,

if : σ = 0 ⇒ σ = σ T (2.49)
Nevertheless, it is argued that under magnetic fields it might be possible to find an asymmetric stress tensor. This could happen because of the way the stress field depends on the electromagnetic field.

Given the Lorentz force on a distribution of charges [START_REF] Espinosa | The magnetic stress tensor in magnetized matter[END_REF] is argued that the asymmetric condition would be given if the microscopic magnetic dipoles in the material are not aligned with the macroscopic field, which will induce a toque per unit volume. Such condition will occur if there is a given interaction not allowing the dipoles to align with the field but in general, it is reasonable to assume aligned dipoles.

#» f = ρ e #» E + #» J × #» B (2.50) Then, #» f = (∇ • #» D) #» E + ∇ × #» H - ∂ #» D ∂t × #» B = (∇ • #» D) #» E - #» B × ∇ × #» H - ∂ #» D ∂t × #» B = (∇ • #» D) #» E - #» D × ∇ × #» E - #» B × ∇ × #» H - ∂ ∂t ( #» D × #» B) = (∇ • #» D) #» E - #» D × ∇ × #» E + (∇ • #» B) #» H - #» B × ∇ × #» H - ∂ ∂t ( #» D × #» B) = ∇ • #» D ⊗ #» E - 1 2 ( #» D • #» E)I + #» B ⊗ #» H - 1 2 ( #» B • #» H)I - ∂ ∂t ( #» D × #» B) = ∇ • σ - ∂ #» S ∂t (

Energy conservation, 1st Thermodynamic Principle:

We consider processes where the energy sources are the mechanical work and input heat. The principle of energy balance establishes that the change rate of the internal energy equals the work done by the volume forces, surface tensions plus the heat that enters the body.

Notation:

Internal Energy density u int Total Energy density

e int = u int + 1 2 #» v • #» v
Heat flux per surface unit #» q Volumetric heat source Q

The rate of change of the total energy is

P tot = P int + P K = D Dt Ω ρu int dΩ + D Dt Ω 1 2 ρ #» v • #» v dΩ (2.52)
The external loads and heat sources /fluxes amount to

P ext + P h = Ω #» v • ρ #» b dΩ + Γ #» v • #» t dS external + Ω ρ Q dΩ - Γ n • #» q dS heat (2.53)
From where

P tot = P ext + P h D Dt Ω ρu int + 1 2 ρ #» v • #» v dΩ = Ω #» v • ρ #» b dΩ + Γ #» v • #» t dS + Ω ρ Q dΩ - Γ n • #» q dS
(2.54) 32

The Solid Mechanics and Heat Transfer Problems

This expression is usually reformulated in three different manners:

Lagrangian ρ D Dt u int = ε : σ + ρ Q -∇ • #» q + #» v • ρ D Dt #» v -ρ #» b -∇ • σ =0 ρ D Dt u int = ε : σ + ρ Q -∇ • #» q (2.55a) Eulerian (a) ρ ∂ ∂t u int + ρ #» v ∇u int = ε : σ + ρ Q -∇ • #» q (2.55b) Eulerian (b) ∂ ∂t (ρe) + ∇ • [ #» v (ρe + p)] = ∇ • (S • #» v ) + ρ #» v • #» b + ρ Q -∇ • #» q (2.55c)

Material physical modelling

An important aspect in the modelling of physical processes is the material characterization that relates the state variables between them. We will first go through the mechanical constitutive relationships that relate the stress state with the strain, strain rate and temperature. We shall present the basic principle of thermal modelling exchange due to energy flux inside a material medium. Finally, we will cover the electromagnetic material description.

Constitutive models

One of the most common ways of describing the mechanical behaviour of metallic materials (not only) is by means of the linear superposition of strain. The basic idea is, as sketched in fig. 2.4 that the total deformation felt by a material volume element is given by the superposition of individual factors. The total stress, on the other hand is the same regardless of the different deformation mechanisms that can be present.

The FORGE modelling tool uses a formulation based on strain-rates, rather than strain (eq. (2.56)). The resolution of the problem is then done in terms of the evolution of the velocity [Mocellin 1999]; The displacements are obtained by time integration.

             ε = εel + εvp + εth εel = D -1 σ εvp = 3 2 # » ε #» σ S , tr εvp = 0 εth = α Ṫ I , tr εth = 3α Ṫ (2.56)
Where εel , εvp and εth represent the elastic, viscoplastic and thermal components of the strain rate tensor. D is the stiffness fourth order tensor, σ is an objective time derivative for the stress tensor. α is the thermal expansion coefficient, Ṫ is the temperature rate of change. ε and σ are the von Mises equivalent strain rate and stress respectively, which are given by

ε = 2 3 ε : ε (2.57) σ = 3 2 S : S (2.58)
Elasticity The classical description of elastic behaviour of solids is given by the Hooke's law

ε el = 1 + ν E σ - ν E tr σ I (2.59)
In terms of the time derivatives we have

εel = 1 + ν E σ - ν E
∇ • #» v = - 1 κ ṗ (2.62)
In the case of isotropic elasticity, computing the trace of eq. (2.60) gives

Tr εel = ∇ • #» v = - 3(1 -2ν) E ṗ ⇒ κ = E 3(1 -2ν) (2.63) For incompressible materials ν = 0.5 ⇒ ∇ • #» v = 0.
Plasticity Once the onset of permanent deformation (Yield stress) has been reached, it is needed to model the stress flow rules. This is an extremely rich field of study in which several descriptions can be obtained depending on the conditions under which we study the materials. As we mention in the previous chapter, the characteristic of the EMF process is that deformation occurs under high strain rates, in the order of 10 3 s -1 to 10 4 s -1 . This requires using flow rules including the effect of strain rate. In this section we just intend to summarize a few of the most used descriptions to model plastic flow under high strain rate conditions.

Johnson-Cook model [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF]]

σ y = A + B εn p 1 + C log ε ε0 1 - T -T 0 T m -T 0 m (2.64)
Where σ y represents the evolution of the yield stress. εp is the equivalent plastic strain, ε the strain rate and T the temperature. A, B, C, n and m are the empirical material parameters. The ε0 parameter was conceived as the effective plastic strain rate at which the strain-hardening parameters A, B and n are determined; today, it is seen mainly as a numerical regularization term, especially when characterizing for a large range of strain rates values. T 0 represents Chapter 2. Physical & Mathematical Framework a reference value for the temperature at which the quasi-static tests are done (assuming that no appreciable self-heating takes place). And finally T m is the fusion or melting temperature.

This model is so far, one of the most used due to its simplicity and versatility. It will also be the one used in our computations in the following chapters.

Hollomon model This is one of the most well known description for work hardening behaviour of metals. Several authors have used this model for the characterization of aluminium alloys in which the strain-rate effects are considered negligible.

σ y = σ 0 εn p (2.65)
Takatsu model [START_REF] Takatsu | High-speed forming of metal sheets by electromagnetic force[END_REF] This flow curve was introduced by Takatsu in order to include the strain-rates effects on the properties of aluminium alloy JIS A1050.

σ y = K 1 ε n 1 p + K 2 ε n 2 p ln ε ε0 (2.66)
Modified path-dependent Zerilli-Armstrong model [START_REF] Gould | A path-dependent constitutive model for gilding copper[END_REF][START_REF] Forde | Symmetrical Taylor impact studies of copper[END_REF] This is a model developed mainly for the characterization of fcc metals, with a special attention to copper. The main characteristic given by the authors is they intend to treat path-dependency in a more physically realistic manner. It

The Solid Mechanics and Heat Transfer Problems

also enables to include in a simple manner the effect of previous loading. Values for the mechanical characterization of some materials using the previously mentioned laws are given in appendix C.

σ y = σ a + Kε n p + µ T µ 0 σ t (ε p , ε, T )φ( ε, T ) (2.67a) µ T µ 0 = 1 -aT (2.67b) σ t (ε p , ε, T ) = η    1 -θ( ε)(α -1) ε p η + 1 - Σ η 1-α 1 1-α    (2.67c) η( ε, T ) = k ε ε0 S( T 1-aT ) (2.67d) φ( ε, T ) = exp [T (-C 3 + C 4 ln( ε))] (2.67e)

Friction

The friction force appears if contact takes place with a relative displacement in the tangential direction besides the normal component. This force opposes the tangential movement between the bodies due to the shear stresses at contact given by

#» t = τ = σ • n -σ n n (2.68)
In the same way the material behaviour laws are needed to relate stress and strain inside a given body, it is necessary to define the relationship between relative displacement and the stress state at the surface. Several laws exist

Chapter 2. Physical & Mathematical Framework

in the literature for that end, here we mention a few of the most known, but before let us define the sliding velocity as

∆ #» v s = ( #» v -#» v tool ) -[( #» v -#» v tool ) • n] n (2.69)
• Power law or viscoplastic friction: This law is derived from the Norton-Hoff viscous material description. It is defined as

τ = -α f K ∆ #» v s q-1 ∆ #» v s (2.70)
Where α f represents the friction coefficient, q the sensibility to sliding velocity and K the material consistency.

• Tresca: It can be obtained as the limit case for the previous model when q = 0. It links the shear stress τ with the plastic flow rule σ y using the tresca friction coefficient

m g = 1 √ 3 σ y m , 0 ≤ m ≤ 1 (2.71a) ∆ #» v s = 0 , if |τ | < g (2.71b) τ = -g ∆ #» v s ∆ #» v s (2.71c)
• Coulomb: It links the critical shear τ c to the contact pressure p = -σ n by means of the friction coefficient µ

τ c = µp ∆ #» v s ∆ #» v s (2.72)
for a more detailed account on the modelling of the contact condition we advise the reading of Wagoner et al. 2001 chapter 6.

Thermal modelling

The thermal problem is obtained from the conservation of energy as shown in eq. (2.55). There are several hypothesis (or models) used in order to get the equations into its explicit representation in terms of the temperature

d dt u int = C p d dt T
where C p is the heat capacity of the body.

• Heat flux or energy diffusion (-∇ • #» q ): The heat flux within a solid material is normally represented by the Fourier relation

#» q = -k∇T (2.73)
Leading to the diffusion term as

-∇ • #» q = ∇ • (k∇T )
where k is the thermal conductivity of the body.

• Mechanical volumetric heat term Qm = ε : σ : Transformation of the dissipated mechanical work into heat. This is the source term from the mechanical solver.

• External volumetric heat sources Qe : In this term we summarize the heat volumetric sources. In electromagnetic applications we have the joule heating term given by

Qe = #» J • #» E (2.74)
where #» J the electric current density and #» E the electric flux field.

The total heat source is then given by

Q = Qm + Qe = ε : σ + #» J • #» E (2.75)

Electromagnetic constitutive laws

Description of the electromagnetic wave propagation phenomena in either free space or material mediums would be incomplete without the constitutive laws relating the different fields. There are three main aspects that have to be introduced as shown in [START_REF] Castro | Mathematical models and numerical simulation in electromagnetism[END_REF].

Chapter 2. Physical & Mathematical Framework

Electric susceptibility and permittivity In dielectric substances, when electric polarization takes place, the typical representation of the electric flux is given by

#» D = (1 + χ e )ε 0 #» E
where χ e is the electric susceptibility. This factor is condensed in the relative electric permittivity by

#» D = ε 0 ε r #» E , where ε r = 1 + χ e (2.

76)

Magnetic susceptibility and permeability Just as in the case of the electric field, magnetic polarization also takes place in several materials when exposed to external magnetic fields. As a parallel to the previous case, the magnetic flux is related to the magnetic field by

#» B = (1 + χ m )µ 0 #» H
where χ m is the magnetic susceptibility. We thus have

#» B = ε 0 ε r #» H , where µ r = 1 + χ m (2.77)
Electrical conductivity (Ohm's law) Most of the known conductive materials used typically exhibit a linear dependence of the electric current density #» J to the applied electric field

#» E. #» J = σ #» E (2.78)
The proportionality parameter σ is called the electrical conductivity. We may also find in literature reference to its inverse value denoted (ρ) electric resistivity.

Summary

At the beginning of the chapter we presented the conceptual coupling between the physics involved in the electromagnetic forming process. Then, we showed the mathematical procedure leading to the mathematical equations relating the different fields. In fig. 2.5 we summarize the main equations describing the whole phenomena and the dependencies among them. This summary allow us to adequately develop the computational model as follows. The subsequent developments are focused on the electromagnetic aspects and the coupling strategies with the thermo-mechanical solver FORGE R . Comme évoqué précédemment, il existe deux méthodes principales concurrentes pour simuler le problème électromagnétique:

Summary

• les éléments finis pour les corps conducteurs et les éléments de frontière pour l'air.

• ou des éléments finis pour tous les matériaux.

Le choix de cette dernière approche comme stratégie de simulation permet de simplifier la manipulation des équations dans un cadre unique -à savoir, les domaines sont différenciés les uns des autres uniquement à travers la définition des paramètres matériau.

En outre, puisque l'approche produit un système linéaire unique pour la résolution du potentiel magnétique, le calcul parallèle est également simplifié. Nous comparons deux stratégies parallèles pour la simulation. La première est basée sur l'utilisation de la bibliothèque externe de calcul numérique PETSc, et la seconde est basée sur l'utilisation d'un solveur interne entièrement codé en Fortran. Il est constaté que le processus d'assemblage des matrices dans le solveur PETSc est très inefficace par rapport au code interne et constitue la principale contrainte.

Dans le but de relier le processus de simulation numérique avec les conditions expérimentales réelles, nous avons ajouté un modèle décrivant le circuit couplé au processus de formage électromagnétique. Cela a été fait par un schéma type circuit électrique qui, couplé au module EF comme une condition aux limites, simplifie le pilotage de la simulation. Nous avons vu que pour analyser les propriétés de stabilité du schéma numérique, il ne suffit pas de se concentrer uniquement sur le système EF, mais de prendre aussi en compte l'influence du couplage avec le modèle du circuit électrique. Cela nous conduit à une méthode d'analyse dans laquelle le système EF est condensé dans une équation représentative scalaire; nous ajoutons aussi l'équation du circuit ainsi que les dérivés en temps des champs concernés. Nous considérons dans les perspectives l'extension de cette analyse à la simulation électromagnétiquethermomécanique.

Les résultats obtenus à partir de la comparaison avec le cas analytique montrent un bon accord avec les résultats numériques. Cette comparaison a également permis de confirmer l'approche quasi-statique pour la simulation des phénomènes de propagation d'ondes électromagnétiques. Un autre résultat important obtenu par cette analyse est qu'une grande partie des pertes énergétiques rencontrées dans le processus peuvent se produire dans le circuit électrique relié à la bobine.

Chapter 3. Computational Modelling of Electromagnetism

We have so far described the theoretical background behind the equations that rule physics in a macroscopic scale and from the continuum perspective.

In this chapter we introduce some of the mathematical concepts behind the finite element formulation needed to transform the strong formulation of the electromagnetic equations (previous chapter) into a discretized weak formulation, needed for modelling in computational terms.

We will then go deeper into the numerical aspects of the electromagnetism.

We will introduce the coupling of the finite element model with a reduced circuit approach of the electromagnetic forming machine needed to pilot the process. This will be followed by a validation case of the electromagnetic aspects only. The last section of the chapter is dedicated to the parallel computing strategy selected in order to increase the performance of the simulation software.

The coupling with the thermo-mechanics will be the main subject of the next chapter; where a special focus is given to the management of motion within (due to) the fully immerse approach selected to discretize the electromagnetic problem.

Electromagnetic Modelling

The potential formulation

In many aspects of physics and several numerical applications, it is more convenient to transform the set of Maxwell equations eq. (2.30) using the potential formulation. Several options could be found as shown in [START_REF] Biro | On the use of the magnetic vector potential in the finite-element analysis of three-dimensional eddy currents[END_REF]].

Nevertheless, the classical φ -#» A formulation is the most widely used. As shown in the previous chapter, the start point are equations eqs. (2.30a) and (2.30d) from where

#» B = ∇ × #» A (3.1) #» E = -∇φ - ∂ #» A ∂t (3.2)

Electromagnetic Modelling

Where #» A is the magnetic vector potential and φ is the scalar vector potential. Using the electromagnetic constitutive laws, we then have

#» D = ε 0 ε r #» E ⇒ #» D = ε 0 ε r -∇φ - ∂ #» A ∂t (3.3) #» H = 1 µ 0 µ r #» B ⇒ #» H = 1 µ 0 µ r ∇ × #» A (3.4) #» J = σ #» E ⇒ #» J = σ -∇φ - ∂ #» A ∂t (3.5)
Introducing these elements into the Maxwell-Ampere equations eq. (2.30b) and the change conservation eq. ( 2.23) it is then obtained

ε 0 ε r ∂ 2 t #» A + σ∂ t #» A + 1 µ 0 ∇ × 1 µ r ∇ × #» A = -σ∇φ -ε 0 ε r ∂ t ∇φ (3.6a) ∇ • (σ∇φ) + ∇ • σ∂ t #» A = 0 (3.6b)
Which is a four variables (φ, A x , A y , A z ) four equations system instead of six variables for a double vector field formulation.

The gauge condition

The double curl in eq. (3.6a) requires conditions to be specified not only on the magnetic potential but also in its first derivative in order to guarantee uniqueness of the solution. One of the most general conditions is the Lorentz gauge

∇ • #» A + 1 c ∂ t φ = 0 (3.7)
Where c is the speed of light in vacuum. This condition is simplified by assuming 1 c ∂ t φ = 0 from where we obtain the Coulomb gauge

∇ • #» A = 0 (3.8)
One of the advantages of such model is that eq. (3.6b) is also simplified. We assume that the differential operators ∇ and ∂ t are commutable, we also assume the electric conductivity σ as uniform in a small domain (corresponding to a single finite element in what follows) then ∇ • (σ∇φ) = 0 (3.9)

Boundary conditions

The electromagnetic problem is going to be solved in a single domain as shown in fig. 3.1. This domain (Ω) is subdivided into three sub-domains: The coil or inductor (Ω I ), the work piece (Ω p ) and the surrounding air (Ω a ). inp and Γ I out . After [START_REF] Biro | On the use of the magnetic vector potential in the finite-element analysis of three-dimensional eddy currents[END_REF]].

The boundary and initial boundary conditions are given as follows

#» A( #» x , 0) and ∂ ∂t #» A( #» x , 0) = 0, #» x ∈ Ω (3.10a) #» B • n = 0 ⇒ n • ∇ × #» A = 0, on Γ B (3.10b) #» H × n = 0 ⇒ ∇ × #» A × n = 0, on Γ H (3.10c) φ = 0, on Γ I out (3.10d) φ = V (t), on Γ I inp (3.10e) Boundary condition 3.10b can be translated into #» A × n = 0 (see 1 )
where n is the outward normal of the surface to which the condition is imposed. This condition is also a symmetry condition for the field #» B. Let's consider a given unit vector belonging to the surface t.

If #» A × n = 0 ⇒ #» A • t = 0 (3.11)
As it will be seen later, this condition corresponds to a Dirichlet condition. On the other hand, condition 3.10c is the natural or Neumann condition emerg-

1 #» B • n = 0 ⇔ ∇ × #» A • n = 0 ⇔ ∇ • ( #» A × n) + #» A • (∇ × n) = 0 ⇒ ∇ • ( #» A × n) = 0 From where one solution is #» A × n = 0. Meaning that imposing this conditions guarantees #» B • n = 0 3.1. Electromagnetic Modelling n ˆA B n  B 0 0 0 ˆ= ⋅ → = × ∇ ⋅ = ⋅ t A A n n B 0 0 ˆ= × × ∇ = × n A n H t ˆt
Fig.

3.2 -Representation of the magnetic boundary conditions on a given elemental face.

ing from the equations. In [START_REF] Biro | On the use of the magnetic vector potential in the finite-element analysis of three-dimensional eddy currents[END_REF] it is said that

∇ × #» A × n = 0 ⇐ #» A • n = 0 (3.12)
We do not need to explicitly use or enforce #» A • n = 0. Nevertheless, it can be noticed that if we accept the continuity of the tangential components on the surface to which this condition is applied, then this condition is a symmetry condition for #» A.

Finite elements for electromagnetism

In traditional FE analysis, we commonly find different elements for which the degrees of freedom are defined at the nodes and the internal values are obtained through interpolation of the nodal values. When using such elements definition for electromagnetic fields several inconvenient may arise. In [Jin 2002] was observed that non-physical or spurious solutions could be found in several cases. This is due to the fact that nodal elements enforce continuity on the normal components of vector fields, which is useful when dealing with divergence or gradient operators. In electromagnetism it is necessary to deal with the curl operator that, in opposition to the divergence or gradient, requires the continuity condition to be applied on the tangential components.

In order to solve this issue, a different kind of elements was introduced by [START_REF] Nédélec | A new family of mixed finite elements in R3[END_REF]. Instead of defining the variables in the nodes, the degrees of freedom for the vector fields are assigned to the edges of the elements.

Nedelec Elements

Lets consider the 2D triangular element given in fig. 3.3. A scalar field defined on the nodes of the element is be given by

φ(t, #» x ) = n φ n (t)ϕ n ( #» x ) (3.13)
For the reference triangle, if we use a P1 (linear) interpolation, the base functions ϕ n ( #» x ) are defined as follows

ϕ 1 ( #» x ) = 1 -x -y ; ϕ 2 ( #» x ) = x ; ϕ 3 ( #» x ) = y
A vector field represented on the edges is defined by

#» A(t, #» x ) = d a d (t)Ψ d ( #» x ) (3.14)
Where a d = #» A• t represents the mean value of the integral of #» A along the edge d, and t is a normal vector tangent to the edge. The vector base functions Ψ d ( #» x ) are computed from the nodal base functions following the convention given in [Jin 2002]

Ψ d ( #» x ) = ϕ i ( #» x )∇ϕ j ( #» x ) -ϕ j ( #» x )∇ϕ i ( #» x ) (3.15)
For the 2D element the vector base interpolation functions are given below and sketched in fig. 3.4. It was explained before that a gauge condition has to be imposed on the system of equations in order to guarantee uniqueness of the solution. One of the advantages of the Nedelec elements formulation is that the Coulomb gauge eq. (3.8) is implicitly guaranteed.

Ψ 1 ( #» x ) = ϕ 1 ( #» x )∇ϕ 2 ( #» x ) -ϕ 2 ( #» x )∇ϕ 1 ( #» x ) = 1 -y, x T Ψ 2 ( #» x ) = ϕ 2 ( #» x )∇ϕ 3 ( #» x ) -ϕ 3 ( #» x )∇ϕ 2 ( #» x ) = -y, x T Ψ 3 ( #» x ) = ϕ 3 ( #» x )∇ϕ 1 ( #» x ) -ϕ 1 ( #» x )∇ϕ 3 ( #» x ) = -y, -1 + x T
∇ • #» A = ∇ • d a d (t)Ψ d ( #» x ) = d a d (t)∇ • Ψ d ( #» x ) ∇ • Ψ d ( #» x ) = [∇ϕ i • ∇ϕ j -∇ϕ j • ∇ϕ i ] + ϕ i ∇ 2 ϕ j -ϕ j ∇ 2 ϕ i (3.16)
The first part of eq. (3.16) is automatically equal to zero and the second part is guaranteed to be null for P1 linear elements. From where

∇ • Ψ d ( #» x ) = 0
It will also be needed to know how to compute the curl of

Ψ d ∇ × Ψ d = 2∇ϕ i × ∇ϕ j (3.17)

Variational approach

Sobolev Spaces The Sobolev spaces are the mathematical vector spaces upon which holds the underlying basis of mathematical analysis and numerical methods for the solution of partial differential equations. These spaces are characterized by allowing the existence of a function (vector or scalar) and its derivative up to a given order. The space also defines a norm for the function and the derivatives. Here, we will mention some of the spaces currently found in finite elements literature.

Let ϕ ∈ R and Ψ ∈ R 3 be a scalar and a vector functions.

Chapter 3. Computational Modelling of Electromagnetism Definition 2. L 2 (Ω) is the space of square integrable functions defined on Ω as

L 2 (Ω) = H 0 (Ω) = ϕ ∈ Ω / Ω |ϕ| 2 < ∞ (3.18)
Definition 3. The set of Sobolev spaces for scalar functions derivable until a degree m H m (Ω) in the domain Ω is given by

H m (Ω) = ϕ ∈ L 2 (Ω) / ∂ α ϕ ∈ L 2 (Ω) ; |α| ≤ m (3.19)
Definition 4. Space of functions vanishing at the boundary

H 1 0 (Ω) ⊂ H 1 (Ω) H 1 0 (Ω) = ϕ ∈ H 1 (Ω) / ϕ = 0 ∈ ∂Ω (3.20)
Definition 5. Space of vector functions with square-integrable divergence

H div (Ω) = Ψ ∈ L 2 (Ω) 3 / ∇ • Ψ ∈ L 2 (Ω) (3.21)
Definition 6. Space of vector functions with square-integrable curl

H curl (Ω) = Ψ ∈ L 2 (Ω) 3 / ∇ × Ψ ∈ L 2 (Ω) 3 (3.22)
Definition 7. Inner products: The following notation for the inner products of the spaces will allow simplifying the notation for the weak forms

Ω f • g dΩ = f, g (3.23) Γ=∂Ω f • g dΩ = f, g | Γ (3.24)
In case of product of tensor fields the , notation implies a double tensor product.

Given the electromagnetic system of equations in potential form

ε 0 ε r ∂ 2 t #» A + σ∂ t #» A + 1 µ 0 ∇ × 1 µ r ∇ × #» A + σ∇φ + ε 0 ε r ∂ t ∇φ = 0 (3.25a) ∇ • (σ∇φ) = 0 (3.25b)
Where #» A ∈ R 3 and φ ∈ R are the field variables of interest, lets define Ψ ∈ H curl and ϕ ∈ H div . Then, the problem consists on finding

#» A ∈ H curl and φ ∈ H div 3.1. Electromagnetic Modelling such that Ψ, ε 0 ε r ∂ 2 t #» A + σ∂ t #» A + 1 µ 0 ∇ × 1 µ r ∇ × #» A + σ∇φ + ε 0 ε r ∂ t ∇φ = 0 (3.26a) ϕ, ∇ • (σ∇φ) = 0 (3.26b)
for all Ψ ∈ H curl and ϕ ∈ H div . From the first Green's vector theorem

Ψ, ∇ × 1 µ r ∇ × #» A = ∇ × Ψ, 1 µ r ∇ × #» A - Γ Ψ • 1 µ r ∇ × #» A × n dS (3.27)
Then, system eq. (3.26) is reformulated as

ε 0 Ψ, ε r ∂ 2 t #» A + Ψ, σ∂ t #» A + 1 µ 0 ∇ × Ψ, 1 µ r ∇ × #» A - 1 µ 0 Γ Ψ • 1 µ r ∇ × #» A × n dS + Ψ, σ∇φ + ε 0 Ψ, ε r ∂ t ∇φ = 0 (3.28a) ∇ϕ, σ∇φ - Γ ϕ (σn • ∇φ)dS = 0 (3.28b)
Regarding the natural conditions we will use the following assumptions

• The electric potential is defined in terms of the input and output potential within the domain of the coil as specified in eq. (3.10)d and eq. (3.10)e.

• No electric flow is allowed on any surface besides Γ I inp and Γ I out .

⇒ Γ ϕ (σ∇φ • n)dS = 0
• No absorbing boundary conditions are used in this work. Thus, the only condition allowing outward flow of the magnetic field is eq. (3.10c)

⇒ 1 µ 0 Γ Ψ • 1 µ r ∇ × #» A × n dS = 0 Finally, ε 0 Ψ, ε r ∂ 2 t #» A + Ψ, σ∂ t #» A + 1 µ 0 ∇ × Ψ, 1 µ r ∇ × #» A + Ψ, σ∇φ + ε 0 Ψ, ε r ∂ t ∇φ = 0 (3.29a) ∇ϕ, σ∇φ = 0 (3.29b)

Chapter 3. Computational Modelling of Electromagnetism

We now define the approximate field solutions corresponding to the finite element discretization of the domain using eqs. (3.13) and (3.14)

φ(t, #» x ) ≈ φ h (t, #» x ) = n φ n (t)ϕ n ( #» x ) #» A(t, #» x ) ≈ #» A h (t, #» x ) = d a d (t)Ψ d ( #» x ) Where #» A h ∈ H rot (Ω h ) and φ h ∈ H div (Ω h ).
Lets define the following notation

M 00 (χ) = Ω ∇ϕ i • χ∇ϕ j dΩ (3.30) M 11 (χ) = Ω Ψ i • χΨ j dΩ (3.31) M 10 (χ) = Ω Ψ i • χ∇ϕ j dΩ (3.32) M cc (χ) = Ω ∇ × Ψ i • χ -1 ∇ × Ψ j dΩ (3.33)
Where χ is a given material property. Furthermore, lets note {a} and {φ} as the tables containing the edge and nodal values a d and ϕ n for the whole domain.

Then, the variational problem can be stated in the discrete matrix form

ε 0 M 11 (ε r )∂ 2 t {a} + M 11 (σ)∂ t {a} + 1 µ 0 M cc (µ r ) {a} +M 10 (σ) {φ} + ε 0 M 10 (ε r )∂ t {φ} = 0 (3.34a) M 00 (σ) {φ} = 0 (3.34b)
At this point we realize from eq. (3.34b) that the resolution of the electric potential field φ does not depend on the magnetic vector potential. This hypothesis relies on the fact that the electric conductivity is only a function of the temperature. This enables to solve the system in a decoupled fashion:

1. Solve the electric potential on the inductor domain only

       M 00 (σ) {φ} = 0 in Ω h I φ i = V on Γ I in φ i = 0 on Γ I out (3.35)
2. Solve the magnetic vector potential in the whole domain. The electric

Electromagnetic Modelling

potential is used as a source term

       ε 0 M 11 (ε r )∂ 2 t {a} + M 11 (σ)∂ t {a} + 1 µ 0 M cc (µ r ) {a} = -M 10 (σ) {φ} -ε 0 M 10 (ε r )∂ t {φ} in Ω h a i = 0 on Γ B (3.36)

Time integration schemes

In eq. (3.35) it is observed that the resolution of the electric potential does not have an explicit dependence on time. Except for the input boundary condition, the distribution of the field does not depend on time. Hence, the finite element resolution of the system may be done once,2 followed by the use of a scaling factor to track the time evolution of the potential due to due boundary condition. This will be seen in the next section.

For the magnetic vector potential eq. ( 3.36) we need to select an appropriate time discretization method. Most publications done in the EMF process domain use a quasi-static assumption. In this work we attempted to verify the validity of this approximation by comparing with the full model and evaluate the influence of the displacement current term.

First, we establish the time-step-wise global linear system to solve. We consider we know the information at a certain time step t and we want to advance to t + ∆t, then

K t {A} t+∆t = F {φ} t+∆t , {A} t , . . . (3.37)
Three time integration schemes are compared

• QSEI: Quasi-static Euler Implicit (One time-step)

K t = µ 0 M 11 (σ) + ∆tM cc (µ r ) (3.38a) F t = -µ 0 ∆tM 10 (σ){φ} t+∆t + µ 0 M 11 {A} t (3.38b)
Chapter 3. Computational Modelling of Electromagnetism

• QSSC: Quasi-static Second-order-Central differences (two time steps)

K t = µ 0 M 11 (σ) + 2 3 ∆tM cc (µ r ) (3.39a) F t = - 2 3 µ 0 ∆tM 10 (σ){φ} t+∆t + 4 3 µ 0 M 11 (σ){A} t - 1 3 µ 0 M 11 (σ){A} t-∆t
(3.39b)

• DYN: Dynamic Newmark (One time-step)

K t =µ 0 ε 0 ∆tb 4 M 11 (ε r ) + µ 0 ∆tb 1 M 11 (σ) + ∆tM cc (µ r ) (3.40a) F t = -µ 0 ∆tM 10 (σ){φ} t+∆t -µ 0 ε 0 M 10 (ε r ) {φ} t+∆t -{φ} t + µ 0 ε 0 ∆tM 11 (ε r ) b 4 {A} t -b 5 {A} t -b 6 {A} t + µ 0 ∆tM 11 (σ) b 1 {A} t -b 2 {A} t -b 3 {A} t b 1 = γ β∆t , b 2 = 1 - γ β , b 3 = 1 - γ 2β ∆t (3.40b) b 4 = 1 β∆t 2 , b 5 = - 1 β∆t , b 6 = 1 - 1 2β (3.40c) β = 1 4 , γ = 1 2 (3.40d)

Machinery modelling and implementation

One of the key points in the simulation of the process consists in including the machinery parameters. This allows ease in the set-up of the simulation by giving the same information as the one used in the real process. The process can be schematized as an electrical circuit as in fig. 3.5. While the electrical architecture of a EMF machine is quite complex, it is often possible to use a simplified model based on 3 equivalent parameters plus the chosen potential (or energy) level. The energy storage capacity is described by the equivalent capacitance C m . Thermal losses in the circuit and connections are summarized into R m and, L m represents the total inductance due to the spatial distribution of the wires.

The machinery is then connected to the work bench, which consists (for our Focusing on the machinery we find that the evolution of the electric potential on the work bench is given by

Machinery modelling and implementation

V (t) = V 0 + 1 C m t 0 i(τ ) dτ + R m i(t) + L m d dt i(t) (3.41)
Where V 0 denotes the initial potential set on the machine and i(t) the electric current running through the main circuit. The latter is a classical output variable of the process measured with a Rogowski coil. This equation is used to define the input electric boundary condition for the inductor FE model by means of the following scheme Equation (3.41) is discretized in time as follows

V (n+1) = V 0 + 1 C m ∆t 2 (i (n+1) + i (n) ) + Int (n) + R m i (n+1) + L m i (n+1) -i (n) ∆t (3.42) Int (n) = t (n) 0 i(τ ) dτ = ∆t 2 n-1 k=1 i (k+1) + i (k)
Since this is an implicit scheme, it requires a predictor-corrector strategy. Let's specify the use of letter n for a given discrete time t (n) and letter k for a given iteration within the implicit scheme when solving the system from n to n + 1.

• At n = 0 and k = 0 we will predict the electric current by approximating Chapter 3. Computational Modelling of Electromagnetism the first time increment of eq. (3.42)

V (1) = V 0 + 1 C m ∆t 2 i (1) + R m i (1) + L m i (1) ∆t from where ⇒ i (n+1) k=0 = i (1) 0 = V (1) -V 0 ∆t 2Cm + R m + Lm ∆t = -χV 0 ∆t 2Cm + R m + Lm ∆t (3.43) • At n = 1 and k = 0 (linear extrapolation) i (n+1) k=0 = 2i (n) -i (n-1) (3.44) • At n > 1 and k = 0 (square extrapolation) i (n+1) k=0 = 3i (n) -3i (n-1) + i (n-2) (3.45)
• For any n and k > 0 the iterative correction is given by a Newton-Raphson computation. We minimize the square of the error function as given by

Error: Φ = 1 2 V (i (n+1) k ) -V (i (n+1) k-1 ) 2 (3.46)
From where

dΦ di k = V (i (n+1) k ) -V (i (n+1) k-1 ) dV di V (i (n+1) k ) -V (i (n+1) k-1 ) 2 i (n+1) k -i (n+1) k-1 (3.47) i (n+1) k+1 = i (n+1) k -θ Φ dΦ di k = i (n+1) k -θ 1 2 i (n+1) k -i (n+1) k-1 (3.48)
χ and θ are purely numerical relaxation parameters to improve the convergence of the numerical scheme. As will be shown later they allow balancing the competition between the finite element model and the circuit equation in the fixed-point scheme. The resolution strategy is summarized in the algorithm ( 1)

Algorithm 1 Electromagnetics FEM-circuit iterative coupling while t < t f inal do Predict i (n+1) k = 1 N orm = 1 Solve ∇ • (σ∇φ) = 0 (eq. (3.35)) while N orm > tolerance & k < N max do Evaluate V (n) = V (i (n-1) , i (n) , i (n+1) k ) Rescale φ with V (n) Solve #» A n+1 = F ( #» A n , #» A n-1 , ∂ t #» A n , ∂ 2 t #» A n , φ n+1 , . . . , φ n-2 ) , (integration scheme) Compute #» J = -σ∇φ -σ∂ t #» A Compute i (n+1) k = #» J • ndS on Γ I inp Correct i (n+1) k+1 N orm = |1 - i (n+1) k i (n+1) k+1 | k = k + 1 end t = t + ∆t end

Validation of the EM finite elements module

For the development and validation stage of the solver, it is recommended to choose a study case that can be easily analyzed through analytical formulations. Fenton 1996 has proposed the ring expansion experiment as the easiest model to perform plasticity characterization on metal specimens submitted to high speed loading.

Description

It consists in two coaxial copper rings as in fig. 3.6. The inner ring serves as the inductor connected to the machine. It will receive the electric loading and generate the repulsive forces on the external ring, causing the latter to deform. For the subsequent analysis, we will consider the following parameters

) ( 1 t i ) ( 2 t i (a) Electric configuration
V 0 [Volts] 5000 C m [µF] 126 R m [mΩ] 48 Case 1 L m [nH] 0 Case 2 L m [nH]
115

1D Analytical Model

By means of Kirchoff's electric circuit laws, we can represent the scheme of fig. 3.5 using the set of equations.

(R m + R 1 )i 1 + d dt [(L m + L 1 )i 1 + M i 2 ] + Q/C m = 0 (3.49a) R 2 i 2 + d dt [L 2 i 2 + M i 1 ] = 0 (3.49b)
Where Q = C m V 0 is the electric charge stored in the main circuit. Considering that no geometry deformation takes place during the EM computation and, assuming the following change of variables.

i 1 = d dt Q = d dt Q 1 (3.50a) i 2 = d dt Q 2 (3.50b)

Validation of the EM finite elements module

The system is now cast as

L m + L 1 M M L 2 Q1 Q2 + R m + R 1 0 0 R 2 Q1 Q2 + 1/C m 0 0 0 Q 1 Q 2 = 0 0 (3.51)
From eq. (3.49) it is possible to say that once the system has stabilized

d dt [L 2 i 2 + M i 1 ] = 0 leading to i 2 = -M L 2 i 1 .
Such assumption allows further simplification of the equations into an equivalent damped resonator given by

Q1 + 2ζω n Q1 + ω 2 n Q 1 = 0 (3.52)
where

ω n = 1 C m L eq (3.53) L eq = L m + L 1 - M 2 L 2 (3.54) ζ = R eq 2 C m L eq (3.55) R eq = R m + R 1 + M 2 L 2 2 R 2 (3.56)
The resistance of an uniform circular section ring is

R = 2a b 2 ρ (3.57)
Where ρ is the electrical resistivity. The mutual inductance is found using Maxwell's principle of geometrical mean distance and the self inductance using Max Wein's formula of linear current distribution [START_REF] Rosa | Bulletin of the Bureau of standards[END_REF]. Leading to

M = µ a 1 + a 2 2 1 + 3 16 2(a 2 -a 1 ) a 2 + a 1 2 ln 4 a 2 + a 1 a 2 -a 1 -2 - 1 16 2(a 2 -a 1 ) a 2 + a 1 2 (3.58a) L i = µa i 1 + 3 8 b i a i 2 ln 8 a i b i - 7 4 -0.092 b i a i 2 (3.58b)

3D Finite Elements Model

In order to find the analytic solution it was supposed that, as well as the ring, the coil is also a closed loop. This allows performing a semi-axisymmetric sim-Chapter 3. Computational Modelling of Electromagnetism ulation in order to save computational time. Prefix semi comes from the fact that an actual 3D simulation is carried out, except for the fact that only a 15 • portion is actually modeled as shown in fig. 3.7. Use of this model requires adapting the input parameters for eq. (3.41) by performing a re-scaling as

V 0 ← 15 360 V 0 , R m ← 15 360 R m , L m ← 15 360 L m , C m ← 360 15 C m (3.59)
For the validation of a FE simulation it is paramount verifying that the final outputs are as independent of the numerical parameters as possible. The main parameters in this case are the mesh size and the time step. Mesh dependency is studied using the 3 mesh refinements levels described in fig. 3.8 and table 3.3.

Stability Analysis

Stability Analysis

Stability and convergence are of paramount importance in the study of numerical simulations. In this section a stability analysis of the time discretization schemes chosen is presented. We follow the analysis technique proposed in [Bathe 1996;[START_REF] Noh | An explicit time integration scheme for the analysis of wave propagations[END_REF]. The discrete-time system is represented as

   x (n) . . . x    t+∆t = M    x (n) . . . x    t + Lr (3.60)
Where x (n) stands for the n-th time derivative vector of the system unknowns, M is the amplification matrix and Lr assembles external loadings and/or previous time steps. Then, stability and accuracy properties can be studied from the characteristics of matrix M From equations eq. (3.6)a and eq. (3.41) two system variables are identified: the magnetic potential vector #» A and the circuit electric current i. The FE model will be represented as a set of scalar "representative" values, from where #» A -→ a. Then, algorithm 1 can be written as

       ä ȧ a ∂ t i i        t+∆t = M        ä ȧ a ∂ t i i        t + . . . (3.61)
Since the selected algorithm is implicit, we will rewrite eq. (3.60) as

B[x] t+∆t = M[x] t .
Then, the characteristic matrix is given by B -1 M. We notice that the method depends on the inversion of matrix B. This resolution does not present a major burden since it is a 5 × 5 matrix. Also, under the parameters definition studied the determinant is different from zero. Meaning that the matrix is not singular.

Identification of the Matrices

• Electric current. i = #» J • ndS = -σ∇φ -σ∂ t #» A • ndS ∼ = -σSGV (t) -σS ȧ from eq. (3.41) ⇒ i (t+∆t) ∼ = -σSG ∆t 2C m + R m + L m ∆t i (t+∆t) + ∆t C m - L m ∆t i (t) + terms (t-∆t) -σS ȧ(t+∆t) (3.62)
Where S amounts for the input surface of the coil and G is a scalar value that accounts for the gradient of the electric potential by -1/l c , l c is a characteristic length.

• Reduction of the FEM system. Lets take the dynamic case as an example.

Stability Analysis

We start by formulating Equation eq. (3.36) in a single element like the one in fig. E.1. The left hand side of the simplified scalar equation is given by

l.h.s = ε 0 ε r m 11 ä + σm 11 ȧ + m cc µ 0 µ r a r.h.s = -σm 10 GV (t) -ε 0 ε r m 10 G V (t) (3.63)
The right hand side is found by first using eq. (3.62) in eq. (3.42), the resulting definition of V (t) is then introduced into the r.h.s of eq. (3.36).

m 11 and m cc are taken as the maximum eigenvalue of eq. (E.7) and (E.9) respectively, and m 10 is linked to chapter E by m 10 = min(M 10 )h (see chapter E for details of the matrices).

The other equations are given by the choice of the time integration method.

1. ä(t+∆t) = b 4 (a (t+∆t) -a (t) ) + b 5 ȧ(t) + b 6 ä(t) 2. ȧ(t+∆t) = b 1 (a (t+∆t) -a (t) ) + b 2 ȧ(t) + b 3 ä(t)
3. FEM simplification (equation eq. (3.63))

4.

∂ t i (t+∆t) = 1 ∆t (i (t+∆t) -i (t) )
5. Electric current (equation eq. (3.62))

Then, all terms related to (t + ∆t) are put on the l.h.s and those related to (t) to the r.h.s. All the coefficients multiplying (ä, ȧ, a, ∂ t i, i) are given in the matrices listed in chapter E.

Results

It is necessary to guarantee that the solution obtained from the numerical model is independent from the numerical parameters. The main parameters to analyse are (a) the relaxation parameters χ and θ, (b) the time-step dt and (c) the mesh size. In this sense, global scalar variables are used to identify and characterize the influence of the numerical parameters over the global simulation.

Chapter 3. Computational Modelling of Electromagnetism

Relaxation parameters

As shown in table 3.1 we have defined two study cases: Case 1 corresponds to L m = 0, which could be a good approximation when dealing with multiple turn coils and, Case 2 L m = 115nH.

The first increment of the electromagnetic simulation is the one requiring more iterations in order to find the dynamic equilibrium state. We evaluate the influence of the relaxation parameters χ and θ during this first increment by evaluating the number of iterations required to reach a relative error in the electric current output of < 10 -6 . From fig. 3.9, it is observed that χ has a quite small influence on the convergence behavior of the algorithm. On the other hand, θ drastically modifies the speed of convergence. It is curious to notice that the curve of convergence changes depending on the value of L m . This actually shows that there exists a competition between the FE system and the circuit model. A small θ means a small weight for the circuit and bigger for the FE output. In fig. 3.9a, L m = 0 so the global impedance of the whole system is basically dominated by the characteristics of the FE model. For such scenario, convergence is improved when θ is small. When L m is different from zero, the competition is more evident as shown in fig. 3.9b. The best convergence is found for θ = 0.5. These results may change from case to case but can be used as a guideline. From here on, we set χ = θ = 0.5.

Stability Analysis

Time step

It is found that, globally, the electric current output follows the behaviour given by the analytic solution. Nevertheless, certain discrepancies are found near the region of the peak current. In fig. 3.10 and fig. 3.11 a zoom around this region is shown.

Chapter 3. Computational Modelling of Electromagnetism

We first remark the impact of the chosen time discretization schemes. While all methods are fully implicit, different degrees of convergence towards the analytic solution are found. As expected, a decrease of the time step generally leads to a better approximation, table 3.4 shows the deviation for the peak current from the analytic solution for each model and each time step.

L = 0nH L = 115nH dt
QSEI QSSC DYN QSEI QSSC DYN 0.1 2.27% 1.14% 0.02% 1.20% 1.31% 0.73% 0.2 3.63% 0.64% 0.56% 2.40% 1.67% 1.63% 0.5 6.89% 1.74% 1.67% 4.46% 3.67% 6.37%

Globally, QSSC and DYN converge faster towards the analytic electric current output solution when compared to QSEI. But, for small time step (0.1µs) spurious oscillations appear in the DYN model when L m = 115nH. These oscillations are amplified when looking at the electric potential in figure 3.13. The period of the oscillations is a regular one corresponding to 0.2µs. Considering the 30cm box surrounding the rings, an EM wave would have bounced back from the walls almost 200 times in this period. The spurious solutions are most likely related to the remnant energy in the waves bouncing from the walls. The imposed boundary condition eq. (3.10c) ( #» H × n = 0) means that we impose the field to be parallel to the outward normal of the walls. Thus, the energy stored in any wave in a direction different from 90 • will be bounced back to the interior of the domain.

Another way to understand the cause of these results is by means of the stability analysis technique proposed in section 3.4. Using the same physical parameters as in the example case the numerical parameters dt = 0.2µs, h = 0.1mm and assuming perfect tetrahedral elements (α = 1) it is found: In order to have a time-integration scheme that is stable, the first requirement is that the norm of the system eigenvalues should be less than one. This guarantees that for null input the solution will always decay to zero. A spectral radii (maximum eigenvalue norm) inferior to one allows filtering the response of higher modes that are spurious responses from the numerical spatial discretization [START_REF] Noh | An explicit time integration scheme for the analysis of wave propagations[END_REF][START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF].

In tables 3.5 and 3.6 we present the eigenvalues resulting from the stability analysis. The first set of eigenvalues are associated to the time-integration scheme of the FE model whilst the second one corresponds to the electric circuit. We first notice that for the Newmark method, the spectral radii is equal to one. It means that any numerical perturbation may (and will) have repercussions on the output response. For instance, the bouncing waves inside the air box.

The second row in tables 3.5 and 3.6 tell us that, for instance, the QSEI approach has more numerical filtering or damping that the other two options. This is the reason why we observed a faster convergence towards the analytic solution for QSSC and DYN as the time-step is reduced.

Looking at the eigenvalues resulting from the circuit parameters, it is evident that the machine inductance has an important weight. For L m = 0nH the eigenvalue (table 3.5) is almost null (0.034) while for L m = 115nH is 0.92 (table 3.6). This is why, when looking at figure 3.12 no spurious oscillations appear, in contrast to 3.13 where the spurious solutions of the dynamic model are highly amplified.

Energy Analysis

It is known that regardless all the advantages of EMF process, it remains quite an inefficient method in terms of energy because the ratio between the input energy (stored in the capacitors) and the actual energy used to achieve the plastic deformation is quite low. [START_REF] Psyk | Electromagnetic forming -A review[END_REF] report that the efficiency rate is between 2% and 25%, the largest is achieved mainly in free expansion cases and the normal expectation should be around 5%. Hence, the importance to characterize the energy distribution and losses throughout the process.

Returning to fig. 3.5 we can establish the following energy balance Chapter 3. Computational Modelling of Electromagnetism

Stability Analysis

E 0 = E Cm + E Rm losses + E Lm Reactive +E available (3.64) 1 2 C m V 2 0 = 1 2 C m V c (t) 2 + R m t 0 i(τ ) 2 dτ + 1 2 L m i(t) 2 + t 0 V (τ )i(τ ) dτ (3.65)
The available energy is then subdivided into

E available = E L 1 + E L 2 + E M Actual usable energy +E R 1 + E R 2 (3.66)

Energy balance at the inductor -work piece level

By means of the analytic model we can understand how the energy is balanced at the inductor -work piece level, we see in fig. 3.14 and fig. 3.15 that the energy is mainly stored as magnetic energy at the beginning of the process, from where the Lorentz forces will be generated. Thermal losses will appear eventually but they become important after the main impulse is given. This energy will contribute to the adiabatic heating of the piece during the deformation process.

Energy balance at the machine level

Using the input parameters for the machinery plus the output electric current and potential from the FE model it is a straightforward computation to obtain the energy balance evolution at the machinery level. Figures 3.16 and 3.17 show the energy balance for case 1 and 2 respectively. It is found that the FE model allows predicting almost the same energy distribution as in the analytic model. This ensures that no numerical diffusion was introduced due to the space or time discretization. The small difference in the behaviour could be attributed to the volumetric effects taken into account by the finite elements discretization. Namely, the wave propagation within the solids, the skin effect and probably, the reflection of wave due to the boundaries of the FE domain.

Parallel computation

3D simulation of coupled electromagnetic/solid mechanics/heat transfer problems usually lead to large number of increments and can be CPU time intensive. Today's computational technology enables to breakdown the problem by means of parallel machines allowing either to reduce the cpu time or increase the level of complexity of the problems we can deal with while preserving reasonable execution times.

Partitioning paradigm

In the following paragraphs we will describe the methodology behind the parallel resolution of the electromagnetic equations:

Mesh partitioning -SPMD

Data parallelism in FORGE and in Matelec follows the Single Program on Multiple Data strategy (SPMD). This method consists in the following idea Each processor runs a complete version of the code, with a data set containing the triangulation of the surface of the tools and a mesh of a sub-domain associated with it. This involves only few change in the data structure of the sequential code, which needs only to deal Chapter 3. Computational Modelling of Electromagnetism with a new type of boundary: The interface between domains [START_REF] Coupez | Parallel simulation of 3D forming processes[END_REF]. The electromagnetic parallelism is an extension of the strategy adopted in FORGE in which, as shown in fig. 3.18, elements are locally known for each domain. Nodes have their own local numbering for each domain and an additional array is created within each partition in order to recognize that a given interface node j in partition p is the same as node k in partition q. Edges follow the nodal structure, we add as well the information of the edge sign at the interface because it might happen that the edge is oriented in the opposite direction in the neighbouring partition.

Matrix and Vectors Assembly

The solver today includes two approaches:

• A FORTRAN layer for data management and pre/post-processing plus a PETSc (C++) [Balay et al. 2014] layer for the actual resolution step.

• A single FORTRAN layer with an in-house built iterative solver based on the Conjugate-Gradient method.

This strategy was adopted because the PETSc library offers a large variety of solution and pre-conditioning methods that can be easily adopted for different kind of problems.

Also, one of the limitations of parallel computing is that as we increase the number of partitions we also increase the number of interfaces. These interfaces imply a communication step. Since communicating a partitioned data can be more time consuming than performing a single operation in a given cpu, we can reach a "saturation point" where the parallel performance starts to decline (see fig. 3.19) as the number of processors increases. This is brought by a competition between the actual computation time versus the communication time, which under the scope of Strong Scalability 3 will rapidly increase. Among the several operations performed within a FE code there are two that are core for the resolution of the involved linear systems. The first is the vector dot product (a = {x} • {y}) and the second one is the matrix-vector product ({y} = A • {x}).

The PETSc library enables to minimize the communication during the matrixvector product by storing matrix A in such a way that each processor knows 3 In parallel computing there are two classifications for its performance: (1)Strong scalability which measures, for a given problem size, the obtained speed up in processing time as the number of cpus increases. (2)Weak scalability which looks at the difference in processing time for proportional increase in problem size and resources. Namely, if we double the resources we also double the problem size and expect the time to be almost the same.

Chapter 3. Computational Modelling of Electromagnetism the full global rows of A. In this way, the operation y i = A ij x j for a given degree of freedom i doesn't need any communication.

In the single Fortran layer setup the matrix is known by parts because each processor only assembles its local elements. This means that we will have partial values y (p)

i that have to be added up at each iteration (y i = p y

(p) i ).
From the previous paragraphs we then see that in PETSc we have an initial communication and assembly work in order to know globally the rows of matrix A ij and avoid subsequent communications due to matrix-vector products. In our in-house model, the initial assembly is always minimal but we do need to perform the communications for the product.

Iterative Solver

For the resolution of the large linear system originated from the EM equations we use an iterative resolution strategy based on the Pre-conditioned Conjugate Gradient (PCCG) method (using diagonal preconditioning).

The Conjugate Gradient (CG) method is an optimization strategy that allows solving problems in the form Ax = b by minimizing the cost function

Φ(x) = 1 2 x T Ax -x T b
In principle, this technique can be applied (only guaranteed to work) for symmetric positive definite (SPD) matrices. One way of enlarging the range of application and optimizing the convergence properties of the method is by using a preconditioner. Let's denote M as the preconditioner matrix that has SPD properties.The conjugate gradient is then applied to the symmetrized preconditioned system

M -1/2 AM -1/2 M 1/2 x = M -1/2 b (3.67)
Instead of using M -1 A that in general can be non-symmetric. The main advantages of this strategy is that (a) we retrieve a symmetric system, meaning that the CG method converge. (b) The spectrum of the matrix

ρ(M -1 A) = ρ M -1/2 AM -1/2 < ρ(A)
meaning that we will improve the convergence rate of the method. An esti-

Parallel computation

mate of the number of iterations required given a certain tolerance value for the error for the CG method is

n CG ∝ ρ(A) 2 (-log ) (3.68) n P CCG ∝ ρ(M -1 A) 2 (-log ) ⇒ n P CCG < n CG (3.69)
For a demonstration of the above statements we recommend the reading of [START_REF] Dolean | An Introduction to Domain Decomposition Methods: algorithms, theory and parallel implementation[END_REF]. For the actual algorithm we do not need to compute M -1/2 but only to apply M -1 as shown in the algorithm 2.

The preconditioning is done by means of the Jacobi preconditioner. It consists on pre-multiplying by the inverse of the matrix diagonal. If we borrow the MATLAB notation, from where diag(A) returns a vector array containing the diagonal terms of matrix A and diag(diag(A)) creates a matrix full of zeros except for the main diagonal filled with the vector's elements then

M = diag(diag(A)). Algorithm 2 Pre-Conditioned Conjugate Gradient (PCCG) Method Initialization r 0 = b -A • x 0 z 0 = M -1 • r 0 p 0 = z 0 relative residual: ρ 0 = 1 i = 0 while ρ i > tolerance and i < N max do ρ i = r i • z i q i = A • p i α i = ρ i p i •q i x i+1 = x i + α i p i r i+1 = r i -α i q i z i+1 = M -1 • r i+1 ρ i+1 = r i+1 • z i+1 β i+1 = ρ i+1 ρ i p i+1 = z i+1 + β i+1 p i i = i + 1 end 3.5. Parallel computation f i = θ i 360 r i r V i = f i V 0 + 1 C m t 0 I(τ )dτ + R m I(t) + L m dI(t) dt (3.70)
Where i represents a given sub-turn of the coil. θ i = 2 • for all. r i is computed from the center of the face of each turn to the central axis and, r = 1 The model is discretized using an isotropic mesh with progressive size refinement. In fig. 3.21 we observe that large elements (approximate size 4mm) are present on the outskirts of the air mesh. The smaller elements are located at the bottom face of the aluminium plate facing the inductor with a size of 0.1mm. This leads to a mesh of 124688 nodes 640940 elements and 792799 edges. The inductors outer borders are meshed with 0.2mm length elements.

N N i r i .
Besides the length size we also set the time step to 0.2µs for this test case. The results shown in fig. 3.22 correspond to the resolution of the first increment of the electromagnetic problem coupled with the circuit model. We focused on this step because it is always the one demanding more cpu time due to assemblies and the fact that since no previous solutions are known the solver takes more iterations between the FE system and the circuit to find the solution. From here we can draw the following results:

• The total time expended to compute the first increment of our simulation is always shorter when using the in-house Fortran solver.

• From the speed up graph we see that for PETSc in 2 and 4 processors, we actually obtain a slow down with respect to a single processor.

• We find an optimum for the in-house code at 16 cores, then we have a slowdown in the performance. On the other hand, PETSc seems to have a more stable growth curve. It would be interesting to test with even more processors.

• The PETSc assembly is extremely expensive. With a single processor we already see that the assembly time is 5times slower than with the inhouse code. 2 procs ⇒ 1300times slower. 32 procs ⇒ 300times slower.

• If we subtract the assembly time from the total computation time we observe that PETSc is actually more efficient up to 12 processors. We would expect the computation to be actually more efficient due to the reduction in cpu time by the lack of communication during the matrixvector product. Nevertheless, the solver time of the in-house code seems to be better than the PETSc code for more than 12 procs.

Once the first increment have been solved, the performance of both codes are comparable. But since in the following we are going to be interested in performing continuous coupling between the electromagnetic resolution and the thermo-mechanical resolution, it means that we will constantly need to pass through the assembly stage of the EM code. We have seen that the assembly is actually the weak aspect of the PETSc interface. From these results, we will favour the use of our in-house code for any application that can be simulated with less than 32 processor.

Conclusions

In this chapter we have reviewed and established the numerical tools needed to model the electromagnetic problem in terms of finite elements analysis. It was observed that due to the presence of the curl operator in the EM equations it is needed to adapt the elements formulation, specifically by using the Nedelec edge elements [START_REF] Nédélec | A new family of mixed finite elements in R3[END_REF]]. They eliminate spurious solutions that would arise if classical Galerkin elements would be used [Jin 2002] and also simplify the definition of the electromagnetic boundary conditions.

Chapter 3. Computational Modelling of Electromagnetism

At the beginning it was said that there exist two main competing methods for simulating the electromagnetic problem: Finite elements for the solid bodies and boundary elements for the air or full finite elements for all materials. The choice of using the latter as our simulation strategy was driven by the simplification of handling the equations within a single framework. Namely, the bodies are different from one another solely by the material parameter definitions. Also, since the approach produces a single linear system for the resolution of the magnetic potential the parallelism is also simplified. We have compared two parallel strategies for the simulation. The first using the well known PETSc library as an external package and the second one with the in-house solver fully coded in Fortran. One of our key findings was that the assembly process of matrices in the PETSc solver was highly inefficient when compared to the in-house code, at the point of putting most of the computational strain in the assembly stage. Our results are limited to the use of 32 processors. It would be necessary to test larger problems in order to evaluate the performance for massive computations

We realized that in order to relate the numerical simulation process to the actual experimental settings it was necessary to include a model describing the machinery behind the electromagnetic forming process. This was done by an electric circuit scheme which coupled with the boundary conditions of the FE module simplifies the piloting of the simulation. Then, we recognized as well that in order to analyse the stability properties of our numerical scheme it was not enough to focus only on the FE system, but also to take into account the influence of the coupling with the circuit model. This led to an analysis scheme where the FE system was compacted into a scalar representative equation, and by input of the circuit equation and also the time derivatives it was possible to study the coupled system. We see as a perspective to enlarge this analysis to the thermo-mechanical/electromagnetic simulation.

The results obtained from the comparison with the analytic case showed a good agreement with the numerical results. It also allowed to confirm the quasi-static approach for the simulation of the electromagnetic wave propagation phenomena. This case also pointed out that a great part of the energetic losses experienced in the process can happen in the electric circuit connected to the coil.

Résumé: Modélisation numérique du problème couplé thermomécanique-électromagnétisme

Dans ce chapitre, nous explorons la méthodologie globale nécessaire pour mettre en place une simulation couplée entre le module électromagnétique et le module thermomécanique. Nous travaillons avec une stratégie de couplage séquentiel qui permet de traiter les différentes physiques de façon individuelle, tout en permettant leur interaction à chaque pas de temps.

L'approche des objets complètement immergés utilisée pour le maillage du problème électromagnétique nécessite des techniques d'adaptation de maillage. Ces dernières permettent de suivre le déplacement et la déformation provenant du calcul mécanique et de l'imposer au maillage de l'air. La première technique implémentée consiste à générer le maillage de façon automatique afin d'adapter la taille des éléments selon les matériaux utilisés. Elle combine le phénomène de la profondeur de peau avec l'utilisation des boîtes de taille de maille qui augmente progressivement. L'objectif principal de cette approche est de concentrer la plupart des éléments autour des solides. Nous avons ajouté le paramètre "distance à la bobine" qui permet de réduire la zone affectée par la taille de maille la plus fine, mais en fonction de la valeur attribuée à ce paramètre la zone affectée est perdue lorsque la pièce s'éloigne de la bobine. Nous avons constaté que cela permet de réduire le nombre des inconnues du système.

Une question à laquelle nous n'avons pas encore répondu est l'effet que cela pourrait avoir sur les résultats physiques et la précision du calcul. Néanmoins, étant donné que l'effet de l'induction diminue rapidement avec la distance, nous ne prévoyons pas un grand impact. Une fois que la force induite a diminué, le principal moteur du processus est l'inertie acquise. La deuxième technique est basée sur la méthode de relocalisation nodale sur la base d'une analogie avec des ressorts non-linéaires. L'un des avantages principaux est qu'il permet de maintenir un maillage unique à travers plusieurs incréments de temps avec un coût minimum. Le temps consacré à résoudre le déplacement des noeuds dans l'air est inférieur à 1% par rapport à un incrément donné.

Une fois le problème de gestion de maillage réglé, reste à résoudre le transport des données en raison du changement topologique dans la discrétisation Dans l'aspect physique du processus, nous avons vu que la simulation par éléments finis a permis de prendre en compte l'effet de la distribution tridimensionnelle du champ magnétique. En fonction du matériau et / ou de la géométrie, cette distribution peut aller dans le sens contraire du processus. La présence des forces de Lorentz sur la face opposée de la pièce de travail peut impliquer un ralentissement de l'accélération de la pièce.

Il est dit dans la bibliographie que le chargement rapide est au centre des avantages que nous pouvons en tirer du processus. Il y a deux phénomènes qui jouent en cette faveur. D'un côté, nous avons le chargement rapide qui se transforme en énergie inertielle. L'inertie contribue au ralentissement et à la redistribution des points de nucléation des fractures au niveau microscopique comme montré par [START_REF] Molinari | Dynamic failure of ductile materials[END_REF]. Nous avons regardé de plus près le second mécanisme qui est lié à l'arrêt rapide contraint du matériau. La dynamique rapide en cas de collision augmente l'effet de l'inertie et les contraintes (moules, matrices) aident à contenir la déformation du matériau.

Ces phénomènes vont de pair avec la réduction du retour élastique, qui est produit par la mécanique d'ondes de choc. Nous pourrions postuler que, si les conditions d'impact sont appropriées pour générer au moins un régime à double vague comme expliqué dans chapter D, le retour élastique pourrait être complètement éliminé.

Chapter 4. Computational Modelling of the Coupled Problem -Objects motion in the immerse FE approach

Increase in the complexity of numerical simulations for forming applications with large amounts of deformation demands inclusion of remeshing or mesh adaptation techniques that allows the fulfillment of the simulation. Extensive documentation exists regarding all the different aspects related to remeshing, mesh adaptation and/ or recovery of fields for mechanical problems, see for instance [START_REF] Zienkiewicz | The finite element method[END_REF]. However, little work has been done regarding the application of these techniques to electromagnetic computations. We found some references in the field of error estimation and mesh adaptation in the works of Zhelezina 2005;[START_REF] Tang | An improved adaptive finite element method for the simulation of electromagnetic field in PPFDs[END_REF]. Most of these works however, deal with static domains and focus on improving the finite elements mesh within a single time increment.

In applications such as electromagnetic forming or induction heating with moving coils, in which electromagnetism is coupled with thermo-mechanics, it is necessary to devise a methodology that enables to follow the movement and/or deformation of the solid parts. One approach relays on the use of boundary elements (BE) for the wave propagation phenomena in the air as in L 'Eplattenier et al. 2009;[START_REF] L'eplattenier | Update on the Electromagnetism Module in LS-DYNA". 12th LS-DYNA Users Conference[END_REF].This eliminates the need for adaptation of the air mesh but increases the cpu cost due to the resolution of the linear system due to the full-matrix nature of BE methods. Also, this approach requires a complex strategy to allow parallel computation as shown in [L' Eplattenier et al. 2010]. Another approach consists on using full finite elements discretization for the whole domain, allowing a simpler definition of the domain in computational terms. In this case it is needed to include a methodology allowing to track the displacement/deformation of the solids that are immersed in the air.

Three main types of mesh adaptation methods are usually described in literature:

• r-adaptation: The mesh topology is preserved and the nodes are relocated in order to smooth the mesh.

• h-adaptation: The mesh is modified by reducing or enlarging the elements size (h) locally.

• p-adaptation: The mesh topology and grid are preserved and the degree of interpolation (p) is modified locally. This technique is mainly used for increasing accuracy but does not tackle the problem of movement,

Solid Mechanics Modelling

meaning therefore it is not interesting in our study.

For the development of our model a combination of h-and r-adaptation are used. As it will be shown, r-adaptation is rather fast with respect to the resolution of the electromagnetic problem (almost negligible) but is limited by the amount of distortion the mesh can endure. h-adaptation allows moving forward at the cost of regenerating the mesh, which induces computational costs related to mesh management and matrix allocations.

In the first part of this chapter we will do a review of the theoretical background of the mechanics and the heat transfer problem, describing how they are handle within the framework of FORGE following the works of [Mocellin 1999;[START_REF] Wagoner | Metal forming analysis[END_REF][START_REF] Pelissou | Discrétisation spatio-temporelle du probleme thermique a deux champs. Application au procédé de forgeage a chaud[END_REF]]. We then discuss the global algorithm that rules the coupling between the electromagnetic module and the thermo-mechanical one. The rest of the chapter is dedicated to the meshing, the nodal relocation strategy and data transport for edge elements, necessary after a remeshing step.

Solid Mechanics Modelling

One way of categorizing solid mechanics simulations is to consider: (a) Structural mechanics resistance, where the main focus is to stay in small deformation regimes. In this case, it is advisable to formulate the continuum and discrete problem in terms of displacements. (b) Forming processes, usually subjected to large deformations of the material for which a formalism in terms of velocities is more adapted. In this section we present the mixed velocity/pressure formulation that improves the stability and convergence requirements in numerical terms [Coupez 1991;Hachem 2009].

We recall that the principal equations defining the flow of matter are given by the conservation of momentum and the conservation of mass

ρ d #» v dt = ∇ • S -∇p + ρ #» b ∇ • #» v = - 1 ρ dρ dt
Where ρ represents the mass density, #» v the velocity field, S is the deviatoric stress tensor and p the hydrostatic pressure, #» b represents the volumetric body As these equations are mathematically and physically incomplete it is needed to include the behavior laws relating the stress and strain states as shown in section 2.4.5, as well as initial and boundary conditions. • On the free surface Γ f r the normal stress vector equals zero: σ • n = 0.

• We can explicitly define the normal stress in Γ t by imposing σ • n = #» t .

• A velocity #» v imp can be imposed on Γ v in conditions such as perfect sticky contact with a rigid tool.

• More generally, the contact with a tool (rigid or deformable) is defined on Γ c . We review the considerations regarding mechanical contact and friction in the next paragraph.

Contact and friction modelling

Computation of contact is done by means of the following steps

Computation of contact distance

Considering a point M belonging to the deformable part, the contact distance is defined positive if the node is outside the tool, zero for perfect contact and negative if penetration occurs, see fig. 4.2a. It is computed as the projection of the vector joining the points M and P onto the outward normal of the body

δ = # » M P • n (4.1)

First order update

Assuming that at time t the contact distance δ t > 0 it is needed to guarantee that δ t+∆t ≥ 0 (fig. 4.2b). If not, nodal displacements should be corrected. The condition is established under two hypotheses: (a) Small time step, from where

δ t+∆t = # » M t+∆t P t+∆t • nt+∆t δ t+∆t δ t + ∆t dδ dt (4.2) ( * ) dδ dt = d dt # » M P + # » M P • dn dt = #» v t tool -#» v t • nt + # » M t P t • dn t dt (4.3) (b) Flat surface during the time step increment dn t dt = 0 ⇒ δ t+∆t δ t + ∆t #» v t tool -#» v t • nt (4.4) Since δ t+∆t ≥ 0 and #» u t ∆t #» v t h( #» u t ) = #» u t -#» u t tool • n -δ t ≤ 0 (4.5)
Chapter 4. Computational Modelling of the Coupled Problem -Objects motion in the immerse FE approach

Signorini or dual Conditions

Verification of contact condition in the master/slave framework typically follows the Signorini conditions. We have already presented the first condition in eq. ( 4.5). The second condition is given by the pressure at contact: The normal stress σ n = σ • n • n is negative if there is contact or equal to zero for a free surface. Hence σ n ≤ 0. The last one -or the dual condition -comes from the fact that either the surface is free σ n = 0 or there is contact h( #» u t ) = 0. So we can summarize the three conditions as

       σ n ≤ 0 on Γ c h( #» u t ) ≤ 0 h( #» u t )σ n = 0 (4.6)

Weak or variational formulation

For the mixed velocity/pressure formulation lets note V the space of admissible velocities, V 0 the space of admissible velocities at zero1 , and P the space of admissible pressures, then

       V = v ∈ (H 1 (Ω)) 3 / (v -v tool ) • n ≤ 0 on Γ c V 0 = {v ∈ V / v • n ≤ 0 on Γ c } P = L 2 (Ω) / Ω pdΩ = 0 (4.7)
We define by means of the principle of virtual work, the variational or weak formulation of the conservation equations

∀(v * , p * ) ∈ V 0 × P    ρ d #» v dt , v * + S, ε(v * ) -p, ∇ • v * - #» t , v * | Γc = ρ #» b , v * ∀v * ∈ V 0 p * , ∇ • #» v + p * , 1 ρ dρ dt = 0 ∀p * ∈ P (4.8)
Depending on the rheology model selected this formulation can be further modified. We will leave it at this stage as the reference equation for the subsequent finite element discretization.

Spatial discretization

At this point we have defined the continuous variational problem. In order to define the discrete form for the finite element formulation we define the subspace Ω h ⊂ Ω which consists in a set of finite elements as shown in fig. 4.3 (h represents the mesh size). The solution of the variational problem is now sought in the approximation subspaces (V h , P h ) where the approximate solution (v h , p h ) is to be found.

The strong coupling between the fields v and p has as consequence that the selection of the approximation subspaces can not be chosen freely, neither independently. In order to guarantee existence and uniqueness of the pair (v h , p h ) it is mandatory to verify the compatibility condition of Brezzi-Babuska.

P + 1 /P 1 elements, Velocity/pressure modeling

The mixed P + 1 /P 1 element (fig. 4.4) was devised as a natural way of satisfying the Brezzi-Babuska condition without the need for artificial stabilization [START_REF] Arnold | A stable finite element for the Stokes equations[END_REF]Coupez 1991]. /P 1 element. The velocity is defined for all the nodes plus the central node (bubble) within the element (P + 1 ) and the pressure is defined only for the nodes (P 1 ).

Chapter 4. Computational Modelling of the Coupled Problem -Objects motion in the immerse FE approach

As for a classical P 1 element, the velocity and pressure fields are interpolated linearly between the nodes of the element. The special characteristic is the addition of an internal node or bubble (as is some times called) to enrich the velocity field. The interpolation function of the bubble is a linear piece-wise function equal to 1 on the central node and equals zero over the external faces of the tetrahedron in 3D, it is continuous within the sub-tetrahedrons. This approach is summarized in eq. ( 4.9) where the velocity is written as the superposition of a linear term #» v l plus the term arising from the bubble #» v b , the pressure only has the linear component.

   #» v h = #» v l + #» v b = N n i=1 #» v li ϕ i + N e k=1 #» v bk ϕ b k p h = N n i=1 p i ϕ i (4.9)
N n represents the total number of nodes and N e the number of elements. ϕ i is a nodal scalar linear interpolation function and ϕ b k the scalar interpolation function for the bubble which is non-zero only within a single element.

The introduction of the bubble term in the element requires to modify the variational formulation as follows.

       ρ d dt #» v l , v * l + ρ d dt #» v b , v * l + S( #» v l + #» v b ), ε(v * l ) -p, ∇ • v * l - #» t , v * l | Γc = ρ #» b , v * l ρ d dt #» v l , v * b + ρ d dt #» v b , v * b + S( #» v l + #» v b ), ε(v * b ) -p, ∇ • v * b - #» t , v * b | Γc = ρ #» b , v * b p * , ∇ • #» v l + p * , ∇ • #» v b + p * , 1 κ ṗ = 0 (
4.10) Since the new term disappears on the faces of the elements #» t , v * b | Γc = 0. Another property imposed on the bubble space is that it has to be orthogonal to the nodal spaces. Thus, [START_REF] Fortin | Experiments with several elements for viscous incompressible flows[END_REF][START_REF] Coupez | Grandes transformations et remaillage automatique[END_REF]]. The resulting system is written as a vector of residual equations

ρ d dt #» v b , v * l = ρ d dt #» v l , v * b = 0 and S( #» v b ), ε(v * l ) = S( #» v l ), ε(v * b ) = 0. Additionally, a static condensation technique is used for eliminating the bubble's inertia term ρ d dt #» v b , v * b [
R = R l , R b , R p t such that        R l ( #» v l , #» v b , p) = 0 = ρ d dt #» v l , v * l + S( #» v l ), ε(v * l ) -p, ∇ • v * l - #» t , v * l | Γc -ρ #» b , v * l R b ( #» v l , #» v b , p) = 0 = S( #» v b ), ε(v * b ) -p, ∇ • v * b -ρ #» b , v * b R p ( #» v l , #» v b , p) = 0 = p * , ∇ • #» v l + p * , ∇ • #» v b + p * , 1 κ ṗ (4.11)
As this problem is highly non-linear due to the material behavior laws, the

Solid Mechanics Modelling

contact and friction problems, it has to be solved using an iterative process like the Newton-Raphson algorithm.

Time discretization

To establish the time integration scheme we define the difference between an increment as the process of advancing from time t to t + ∆t and an iteration as the process of correcting a given solution by trial and error. A sequence of iterations will take place at time t + ∆t at which the variables are unknown fig. 4.5. From eq. ( 4.11) the condition R = 0 has to be enforced at any given time. Depending on the time discretization used in eq. ( 4.11) the residual will contain terms from time t + ∆t as well as from previous time steps. Namely, R = R(t + ∆t) + R(t). For instance, the time derivatives of the velocities are approximates as

t t t t ∆ + 1 2 3 i f ) (t f ) ( t t f ∆ +
ρ d #» v l dt = ρ ∆t #» v (t+∆t) l -#» v t l (4.12)
Using a Taylor expansion between a known iteration i and the unknown i + 1 we can write

R (i+1) = R (i) + ∂R ∂y (i+1) δy (i+1) 0 = R (i) (t + ∆t) + R(t) + ∂ (R(t + ∆t) + R(t)) ∂y (i+1) δy (i+1) ∂R(t + ∆t) ∂y (i+1) δy (i+1) = -R (i) (t + ∆t) + R(t) (4.13)
Where y stands for any of the principal variables in the system ( #» v l , #» v b , p). Defin- at each iteration the following matrix system has to be solved

   K ll K lb K lp K bl K bb K bp K pl K pb K pp       δ #» v l δ #» v b δp    = -    R l R b R p    (4.15)
The final solutions are obtained by the progressive update

       #» v (i+1) l = #» v (i) l + δ #» v l #» v (i+1) b = #» v (i) b + δ #» v b p (i+1) = p (i) + δp (4.16)
Several simplifications can be done to system eq. ( 4.15). We already saw that K lb = K bl = 0 due to the orthogonality property. Since the bubble is defined within a single element a local condensation can be done in order to reduce the number of unknowns. Line 2 of system eq. ( 4.15) implies

δ #» v b = -(K bb ) -1 (R b + K bp δp) (4.17)
Reintroducing it into the linear system we obtain the reduced representation

K ll K lp K pl -K pb (K bb ) -1 K bp -K pp δ #» v l δp = - R l R p -K pb (K bb ) -1 R b (4.18)
Note: for incompressible materials K pp = 0. This term is different from zero in the case of compressible elastic deformation in solids. The final system obtained here is then solved by preconditioned iterative solvers using PETSc [Balay et al. 2014].

Heat Transfer Modelling

While the heat transfer phenomena is not a main source for the loading in the electromagnetic forming process, its influence can not be neglected because the short process time does not enable diffusion to take place. Since a

Heat Transfer Modelling

great part of the energy involved in the process is recovered and lost as heat, the temperature accumulated within the work piece will modify its material properties. It may lead to higher mechanical deformation and also modify its electromagnetic properties such as the electric conductivity or magnetic permeability.

In this section we summarize the current modelling strategy for the heat transfer phenomenon. As for the mechanical fields, the resolution consists on a single field formulation, where the temperature is the principal unknown of the system.

The heat transfer problem is characterized by the energy conservation principle

ρC Ṫ + ∇ • #» q = ẇ (4.19)
Where ρ is the material density, C the specific heat, #» q the heat flux, Ṫ is the total time derivative of the temperature in a lagrangian framework. ẇ is the heat source or function of heat generation. The heat flux inside the body is related to the temperature by means of the Fourier conduction law #» q = -K∇T (4.20)

The strong formulation is then written as ρC Ṫ -∇ • (K∇T ) = ẇ (4.21)

Boundary Conditions

Considering the solid domain Ω of fig. 4.6a and its external boundary

Γ = Γ 1 ∪ Γ 2 ∪ Γ 3 ∪ Γ 4 .
The main boundary conditions that can be defined are

• On Γ 1 the temperature is imposed to a known value T = T imp .

• On Γ 2 the heat flux is imposed to a known value -K∇T • n = φ imp .

• Pure convection and radiation are mixed in a single convection coefficient h by saying h = h cv + h r . Where h cv will depend on external conditions of the air or fluid in contact and hr ε r σ r (T + T ext )(T 2 + T 2 ext ). ε r represents the surface emissivity and σ r = 5.67 • Contact with another solid body will imply a resistive surface contact. For this, an equivalent convection coefficient for surface conduction h cd is used. It takes into account the conductivities of the two bodies plus the surface rugosity that traps air in-between.

× 10 -8 [W m -2 K] the Stefan-Boltzmann constant. Then -K∇T • n = h(T -T ext ) motion in the immerse FE approach Ω { } imp n q φ = ⋅ Γ : 2 { } ) ( : 2 3 T T h q cd - = Γ { } ) ( : 4 ext T T h q - = Γ { } imp T T = Γ : 1 (a) Boundaries

Weak formulation

The resolution of the heat transfer problem eq. ( 4.21) is done by projecting into the space V T = H 1 0 of test functions ϕ T for the temperature. Then

ρC Ṫ , ϕ T + K∇T, ∇ϕ T + h cd T, ϕ T | Γ 3 + hT, ϕ T | Γ 4 = ẇ, ϕ T -φ imp , ϕ T | Γ 2 + h cd T 2 , ϕ T | Γ 3 + hT ext , ϕ T | Γ 4 (4.22)

P1 elements

The temperature field is found by discretizing the domain with tetrahedral P 1 elements having the temperature defined on the nodes of the element as in fig. 4.7.

Heat Transfer Modelling

The discrete formulation of the time at any given time t and geometrical point

x is then given by the Galerkin interpolation

T (x, t) = N n i=1 T i (t)ϕ i (x) (4.23)
Where N n represents the number of nodes, ϕ i are the linear interpolation functions for the temperature.

Introducing the interpolation functions into eq. ( 4.22) we can retrieve the following system of equations

M ∂ ∂t {T } + K {T } = {Q} (4.24)
Where {T } represent the set of nodal values of the temperature. The matrices C and K and the source term vector {Q} are defined as The matrix differential system eq. ( 4.24) is non-linear because material coefficients and the surface integral terms depend on the temperature solution. the resolution is done following an implicit scheme for the time integration as explained below

M = N elts e=1 Ωe ρCϕ T j ϕ T i dΩ e ( 4

Time evolution

System eq. ( 4.32) is discretized in time by means of a two-time-steps method as follows:

• Equilibrium is computed at a time t * such that

t * = α 1 t n-1 + α 2 t n + α 3 t n+1
and ∆t 0 = t nt n-1 and ∆t = t n+1t n . The parameters α 1 , α 2 , α 3 are selected such that α 1 + α 2 + α 3 = 1.

Chapter 4. Computational Modelling of the Coupled Problem -Objects motion in the immerse FE approach

• The unknown field is also interpolated as

U * = α 1 U n-1 + α 2 U n + α 3 U n+1 (4.28)
• The time derivative of the field is computed as

∂ ∂t U * = β 1 U n-1 + β 2 U n ∆t 0 + γ 1 U n + γ 2 U n ∆t (4.29) With γ 1 = -γ 2 and β 1 = -β 2 = γ 2 -1.
• Unconditional stability and consistency are guarantee using the relations

γ 2 ≥ 1 2 & α 1 > 1 2 (1 -γ 2 ) (4.30)    α 2 = 1 -α 1 + ( 1 2 -α 1 ) ∆t 0 ∆t -γ 2 2 (1 + ∆t 0 ∆t ) α 3 = (α 1 -1 2 ) ∆t 0 ∆t + γ 2 2 (1 + ∆t 0 ∆t ) (4.31)
• The system matrices and the vector are also described at time t * using the same interpolation as in eq. ( 4.28) so that

M * ∂ ∂t U * + K * U * = F * (4.32)
By substitution of the time derivative and the field description at t * it is found

M * γ 2 α 3 ∆t 2 + K * U * = F * + M * α 1 γ 2 α 3 ∆t - β 1 ∆t 0 U n-1 + α 2 γ 2 α 3 ∆t + β 1 ∆t 0 + γ 2 ∆t U n (4.33)
And the nodal unknowns at time t n+1 are recovered by

U n+1 = 1 α 3 [U * -(α 1 U n-1 + α 2 U n )] (4.34)
This method allows selecting several sets of parameters that will derive in different numerical performances. For instance, the well known Crank-Nicholson method would be found for α 1 = 0, α 2 = 1/2, α 3 = 1/2, γ 2 = 1, β 1 = 0. The method that proved the best quality of results in the work of Pelissou 2005 was the Implicit Dupont which is obtained for the set α 1 = 0, α 2 = 0,

α 3 = 1, γ 2 = 4.3. Global algorithm 3/2, β 1 = 1/2.

Global algorithm

Multi-physical coupling is a challenging process since the different specific space and time lengths require a proper strategy enabling the interaction between the different aspects (heat-transfer, mechanics, electromagnetism). We need to distinguish between two concepts. The first one is related to the resolution of the linear systems arising from each physics. The second one is related to the time evolution of the variables of each system and how we allow them to interact.

• Resolution Strategy

Lets assume we have n systems, each one composed of several state variables {X} i (where i is the system under consideration). Each system represents a specific phenomenon on a given physical system. The coupling between these phenomena can be categorized as

Strong Coupling    [A] 11 • • • [B] 1n . . . [B] 1n • • • [A] nn       {X} 1 . . . {X} n    =    {f } 1 . . . {f } n   
All physics are included in a single solver or framework. The main advantage is the accuracy of the results. Its disadvantage comes in the computational burden to generate the large matrices describing all the physics together. An application in which this kind of coupling can be necessary is the resolution of magneto-hydrodynamics (see for instance [Badia et al. 2012]).

Chapter 4. Computational Modelling of the Coupled Problem -Objects motion in the immerse FE approach

Weak coupling

[A] 11 {X} 1 = {f } 1 - j =1 [B] 1j {X} j . . . [A] nn {X} n = {f } n - j =n [B] 1j {X} j
Each physics is solved individually, coupling is achieved by updating the inputs of the current system with the solutions found from the previously solved systems.

• Time evolution

Direct coupling:

It can be applied to either Strong or weak coupling. All system variables advance together. Several iterations within a single time step will be required for this approach.

Sequential coupling:

Applies only to the weak coupling scheme. The idea is that during a given time step increment each system is solved once. The communication is performed at the end of each resolution but no correction step is required.

Loose coupling:

Similar to sequential coupling. The difference being that communication from the electromagnetic resolution is done at the end of the full time period. This means that for the EM simulation no geometry evolution is taken into account.

We have selected the sequential approach because it allows a clearer and adapted management of the simulation without highly affecting the accuracy [START_REF] Alves | Magnetic Pulse Forming: Simulation and Experiments for high-speed forming processes[END_REF]]. Also, it enables taking advantage of the existing technology within the FORGE software for automatic mesh generation as will be seen after.

Global algorithm

Automatic pre-meshing using the skin effect

Using the automatic mesh generation capabilities included within FORGE makes it is possible to use the material and process characteristics in order to generate adapted meshes. In this case we consider the skin effect. This phenomenon appears when a time-varying current circulates through a solid conductor, generating magnetic fields that induce eddy currents within the core of the conductor expelling the electric flow to the outer surface as shown in fig. 4.10. In a harmonic oscillating regime the skin depth is given by

δ ω = 2ρ ωµ r µ 0 (4.35)
Where δ ω is the depth from the surface in which the current density has dropped to 0.37(1/e) times its value at the surface. ρ is the material electric resistivity.

ω represents the electric signal frequency. µ 0 is the magnetic permeability of void and µ r the relative value for each material. In order to get a time evolution estimate, the following approximation can be done: Consider the first quarter of a pulsed signal so that t ∈ [0; T /4] where T is the signal period.

δ t = 2ρ µ r µ 0 T 2π = ρT µ r 4π10 -7 = ρT 4π 2 10 -7 µ r δ t [mm] = 10 6 2π 10 ρ µ r √ t = 10 6 2π 10 ρ µ r √ N √ ∆t = f m ∆t[s] (4.36)

R-adaptation

The element size in the selected domain is then estimated by h = δ t N elts/δ (4.37)

Where N elts/δ is a pre-defined number of elements in the skin depth. Once all the parameters are set, this algorithm is used to automatically generate an adapted mesh at the beginning of the simulation setting process as well as it can be automatically activated during the simulation to regenerate the mesh as it will be shown.

R-adaptation

One of the advantages of these type of algorithms is that there is no new mesh management associated. The domain discretization is kept fix while the adaptation is introduced by virtual displacements . The drawback is that the applicability of this method is limited by the obvious distortion of the elements. To bypass this problem we verify the mesh at the end of the adaptation. If the motion in the immerse FE approach quality is not verified, a new mesh will be generated for the subsequent timestep increment, which is the object of the next section. Several algorithms can be found in literature to induce virtual displacements in a mesh given the displacements of some of its boundaries. The most classic ones are the Laplacian smoothing [Herrmann 1976] methods or the barycentric relocations. The latter allows choosing how the weights associated to the neighbouring nodes are computed, in case the weights are simply given by 1/N where N represents the number of neighbouring nodes, we get an the iterative form of the Laplacian smoothing. Another type of algorithm is based on virtual springs. In this methods we consider that the elements edges are elastic springs and the problem could be reduced to a structural mechanics problem of elastic beams. Several variants and refinements of this technique can be found, see for instance [START_REF] Farhat | Torsional springs for two-dimensional dynamic unstructured fluid meshes[END_REF][START_REF] Degand | A three-dimensional torsional spring analogy Bibliography method for unstructured dynamic meshes[END_REF][START_REF] Bottasso | The ball-vertex method: a new simple spring analogy method for unstructured dynamic meshes[END_REF][START_REF] Zeng | A semi-torsional spring analogy model for updating unstructured meshes in 3D moving domains[END_REF]]. Lets take a set of springs connected to a single node p as shown in figure 4.13. The force on p due to the given displacement in node q Fpq is given by

Non-linear springs

R-adaptation

Fpq =    k x pq 0 0 0 k y pq 0 0 0 k z pq    (ū p -ūp ) (4.38)
Where k d pq are the non-linear stiffness coefficients defined for each direction so as to enable anisotropic behaviour. ūi represent the displacement of the nodes.

The equilibrium is reached when for all nodes.

q =p F d pq = 0 ⇔ 0 = q =p k d pq u d p -u d q (4.39)
The non-linear stiffness of the springs is defined as in (4.40)

k d,i pq = u d,i p -u d,i q x d,t p -x d,t q (4.40)
This strategy seeks to minimize the average deformation by continuously modifying the stiffness of the edges. When a given element gets excessively distorted it will get extremely rigid and thus force to propagate the deformation towards the neighbouring nodes.

Solution strategy

Given the non-linear nature of the method we need to chose an iterative algorithm in order to solve it. For the fixed-point resolution of the system it was chosen to use a relaxed approach using the following scheme

u d,i+1 p = (1 -θ)u d,i p + θ q =p k d,i pq q =p k d,i pq u d,i q (4.41)
Where θ is the relaxation parameter that we have defined as θ = 2 1+ √

(Inverse of the golden ratio).

It was also selected to work with a matrix-free approach in order to reduce the memory allocation time. For the iterative resolution of the algorithm in parallel computation we simply need to define Another reason to prefer the iterative-matrix-free method was the easiness to include symmetry planes constraints into the solver in any specific direction as described next.

Non-constrained symmetry planes

The use of symmetry planes is of great advantage in many FE analyses since it allows reducing the computational time when the geometry and the physics of the problem allows it. When a symmetry plane matches with one of the main directions X, Y or Z we can easily block the displacement by setting u d = 0. When this is not the case, we can further correct the displacement field by the iterative approach

ūcorrected p = ūi+1 p -ūi+1 p • n n (4.44)
Otherwise, it would be required to assemble the matrices including all the degrees of freedom of the system and include rotation matrices in order to transform the plane constraint into a local system for the element to which the condition is applied.

Field transport in electromagnetism

Lets consider the two meshes given in figure 4.16, both meshes are supposed to have the same boundaries while not the same internal topology. Additionally, lets consider a given vector field which is (should be) identical in both meshes. When the field is represented by point data (on the nodes or integration points) it is straightforward to transport the field from the old triangulation (mesh) T 1 h to the new one T 2 h by means of an interpolation approach or a recovery strategy. The edge data available due to the use of Nédélec elements introduces the difficulty that a single edge in T 2 h may run through several elements in T 1 h . Since it is not a point data, the question is, how should it be interpolated? Or averaged? We shall first define mapping as the whole process of representing the edge data in one mesh as edge data in another mesh. This mapping process is divided in three stages categorized in two types of processes as represented in figure 4.17

Field Transformation -element wise

Once the different stages have been identified, it is possible to define the mathematical expressions for the field transformation. Let #» A( #» x ) represent the field in the original triangulation T 1 h and

#» A * ( #» x ) the field in the new triangulation T 2
h . Following the notation from the previous section in eq. (3.14) let Ψ( #» x ) represent the base interpolation function for the edges.

The problem statement consists in finding

#» A * ( #» x ) #» A( #» x ) (4.45)
Which can also be represented by the minimization problem

min #» A * ( #» x ) - #» A( #» x ) 2 (4.46)
Since the field is known in a discrete form

#» A * ( #» x ) = j a * j Ψ * j ; #» A( #» x ) = i a i Ψ i (4.47)
Where a i = #» A • t is the value of the line integral of #» A over the edge (i). t is the tangential unit vector of the edge. (i) and (j) are indexes counting the edges on T 1 h and T 2 h respectively. Since we are interested in finding

#» A * , the minimization problem is found by projecting #» A * ( #» x ) - #» A( #» x )
on the base functions of T 2 h as: 110 4.6. Field transport in electromagnetism

j a * j Ψ * j - i a i Ψ i , Ψ * k = 0 (4.48) j Ψ * j , Ψ * k a * j = i a i Ψ i , Ψ * k (4.49)
Obtaining the linear system

M {a} * = b({a}) (4.50) where        M jk = Ψ * j , Ψ * k System matrix {a} * j = #» A • t * j unknowns b k ({a}) = i a i Ψ i , Ψ * k Loading (4.51)
Notice that the terms a j Ψ i are known on the triangulation T 1 h while Ψ * k is evaluated on T 2 h . This makes it necessary to include the interpolation or transport stage before assembling the loading in eq. ( 4.51). The weights ω k associated to the integration points of the element in the old mesh depends on the relative coordinates of the point p: (η p , γ p , ξ p ) by

Field Interpolation

ω 1 =c (1 -η p -γ p -ξ p ) + d (η p + γ p + ξ p ) (4.54a
)

ω 2 =cη p + d (1 -η p ) (4.54b
)

ω 3 =cγ p + d (1 -γ p ) (4.54c
)

ω 4 =cξ p + d (1 -ξ p ) (4.54d) a = 1 4 1 - 1 √ 5 (4.55a) b = 1 -3a (4.55b) c = -2a -b 3a 2 -2ab -b 2 (4.55c) d = a 3a 2 -2ab -b 2 (4.55d)
As shown before, this operation is local in the sense that it only requires data from a single element. It has the advantage as well of not requiring any special treatment for the interfaces between bodies because each element belongs to a single body. Lets remember that in the interface between a metallic body and a non-conductor body a discontinuity of the normal component of the electromagnetic field may be present (only tangential continuity is enforced).

Field Transformation -node wise

It could also be proposed to use nodal data for the transport of the field but there are two main inconvenients: (1) Nodal vectors would render difficult handling discontinuities between interfaces and (2) definition of the boundary conditions are also complicated in term of nodal variables. Nevertheless, a scheme going back and forth between nodes and edges was tried out and compared with the previously mentioned method.

The node-to-edge method is described as follows:

Ā(x n ) = N edges j=1 a j Ψ j (x n ) (4.56) a d = 1 N nodes N nodes i=1 t • Ā(x n ) (4.57)
Nodal values Ā(x n ) are computed by adding up the contribution of all the edges converging on the node (eq. ( 4.56)). Edge values a d are approximated by

Field transport in electromagnetism

averaging the projections of the nodal vectors on the edge (eq. ( 4.57)). N nodes is the number of nodes belonging to a single edge. The test case is a simple unit squared 2D mesh as shown in fig. 4.19. Each element contains three basic nodes and three additional virtual nodes in order to verify whether the data diffusion of the mapping can be reduced. In fig. 4.17 this is noted as 1st or 2nd order accuracy mapping. This category is also applied for the edge-element mapping, in which 1st order refers to the use of a single integration point by element while 2nd order refers to the use of the 2nd order integration points.

Validation of the transport methodology

The interest behind this case is that it is possible to easily visualize the operators behind the mapping process. Assuming that some edge data is known, it is desired to know how the application of these methods will affect the edge data.

    . . . a d . . .     new = Ω Ψ i • Ψ j dΩ -1     . . . • • • . . . ω int Ψ x ω int Ψ y ω int Ψ z . . . • • • . . .     OP El2Ed (Transport)     . . . • • • . . . N edges j=1 Ψ j (x int ) . . . • • • . . .     OP Ed2El (Transform)     . . . a d . . .     old (4.58)
The edge-to-element operator (OPEd2El) and the element-to-edge recovery Chapter 4. Computational Modelling of the Coupled Problem -Objects motion in the immerse FE approach operator (OPEl2Ed) are shown in eq. ( 4.58). They represent the processes of transforming the edge data into vectors defined at the integration points (OPEd2El) and then proceed with the transport or interpolation to the new mesh. The weight ω int correspond to the ones defined in eq. ( 4.54). Besides the operators it is also needed to consider the matrix from the linear system upon which the mapping process depends. 4.59 summarizes the transformation stages by means of the edgeto-node operator (OPE2N) and the node-to-edge recovery operator (OPN2E).

    . . . a d . . .     new =     . . . • • • . . . 1 N t d x 1 N t d y 1 N t d z . . . • • • . . .     OP N 2E     . . . • • • . . . N edges j=1 Ψ j (x n ) . . . • • • . . .     OP E2N     . . . a d . . .     old (4.59) Equation
We recall that a projection of the data into original mesh should return an identity matrix operator. From this statement it is interesting to see that for the 1st order scheme of the edge-element mapping noise is induced onto the data, This is evident from the resulting dense matrix. From the diagonal values we see that each element of the data array multiplies its own value, and every element has an influence on the rest. This scheme is obviously not suited for the process.

1st Order 2nd Order

Ed2El        2.2 -0.8 -1.41 0.2 0.2 -0.8 2.2 -1.41 0.2 0.2 -0.35 -0.35 2 0.35 0.35 0.2 0.2 1.41 2.2 -0.8 0.2 0.2 1.41 -0.8 2.2               1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1        (4.60)
On the other hand, use of the 2nd order approach returns a perfect identity matrix, meaning that the transformations did not induce any data diffusion. The difference between both methods lies on the assembly process of the linear system. The solution of the finite element model returns a linear field on the elements (P1) but, the degree of interpolation is reduced to P0 if the data is stored in a single integration point. Thus the noise found after retrieval as a

Data research in parallel computation

P1 data set on the same mesh.

1st Order 2nd Order

E2N        1 0 -0.71 0 0.5 0 1 -0.71 0.5 0 0 0 2 0 0 0 0.5 0.71 1 0 0.5 0.2 0.71 0 1               1 0 -0.47 0 1/3 0 1 -0.47 1/3 0 0 0 2 0 0 0 1/3 0.47 1 0 1/3 0.2 0.47 0 1        (4.61)
In the case of nodes to edges it is observed that in both cases the mapping is different from identity. This means that numerical noise will be introduced in the target field. Looking at the non-diagonal terms it is seen that the increase in the number of nodes did slightly reduce the noise but no remarkable improvement is obtained.

Data research in parallel computation

Notation: In this section, whenever we refer to processor it will be taken as a synonymous for a machine cpus core.

When we proceed with a mesh regeneration algorithm in our current solver there are no means for controlling the location of of the mesh partitioning with respect to the processors.

To clarify the idea, let us look at figure 4.20. In this idealized case we have an initial square domain partitioned in four squared sub-domains belonging to the processors (0,1,2 & 3). At a given moment it is decided to regenerate the mesh. The hypothetical result is the square divided in four triangles. So not only have the processors changed their "geometrical location" with respect to the original domain, but also the sub-domains are modified.

Then, before applying the mapping strategy, it is needed to locate for each new processor the processor from the old partitioning that will contribute to the data transfer. If we take process 3 for the triangle-partitioned configuration in figure 4.20, then, the corresponding old partitions will be p0 and p2.

Domain Hierarchy: Shape resemblance ranking

In order to establish in which order we shall start the search of data between the new partitioning and the old ones we are going to create a classification by means of a resemblance factor.

Coverage Ratio

Given a partitioning domain p in the new increment and the is partitioning from the old increment we define the Coverage Ratio of i over p (CR p i ) as

CR p i =        0 ξ ri -ξ rp ≤ 0 1 2 + 1 4 (ξ ri -ξ rp ) 0 < ξ ri -ξ rp ≤ 2 1 2 < ξ ri -ξ rp

Data research in parallel computation

Where ξ is a reduced coordinate between the origins of the circumscribing spheres of each domain as shown in figure 4.21. ri and rp are subscripts to indicate the zone covered from the center to the outer limit of the circumscribing sphere of domain i or p respectively.

Resemblance Factor

The resemblance factor will be given as a combination of the Coverage Ratio and the distance between the centers of the spheres by

RF p i = 1 2 exp -D 2 ip + exp -(1-CR p i ) 2 (4.63)
The ranking is then performed, for each partition p, by ordering the is domains with respect to RF p i . To illustrate the process lets take the following example: The partitions shown in figure 4.22 are characterized by their circumscribing spheres (see tables 4.1). The complete procedure is defined in algorithm (3). CR RF P0 0.49 0.98 0.89 0.49 0.98 0.89 0.26 1.00 0.97 0.26 1.00 0.97 P1 0.49 0.98 0.89 0.26 1.00 0.97 0.26 1.00 0.97 0.48 0.98 0.89 P2 0.26 1.00 0.97 0.48 0.98 0.89 0.48 0.98 0.89 0.26 1.00 0.97 P3 0.26 1.00 0.97 0.26 1.00 0.97 0.48 0.98 0.89 0.48 0.98 0.89

The results for the 2D case are shown in table 4.2. Then, the ranking process follows simply by comparison of the results of the resemblance factor. If we take for instance p2 from the new partition, then from RF the research follows the following order (p0, p1, p2, p3). Additionally, if all the elements have been found at a given pi it is not necessary to continue the search on the rest of the domains, so the research will focus only on (p0, p1).

Computational Time Analysis

In order to evaluate the impact of each strategy included we are going to compare the following cases: 1. Remesh if there are excessively distorted elements (ratio between the minimum and maximum edge lengths less than a given tolerance).

2. Remesh if there are elements with negative volume. The model used consists on two parallel rectangular solids as shown in fig. 4.11.

   Q e = min(h i ) max(h i ) < tol = 1% ∃ Ω e < 0 ( 4
We assign copper as material to both parts. The FE mesh in composed of 46447 nodes, 254298 elements and 307396 edges. The time step is fixed at 0.2µs. The total time consumed for computing 50 increments is always bigger in the case of forced meshing. We also observe a loss of performance when going from 4 to 8 processors. These results are related to the fact that the remeshing and preparation are required several times for the forced case. This also implies more calls to the transport and interpolation routines.

Remeshing and Preparation

Following the different steps in a given time-increment, the first two stages to intervene are the remeshing (if it was activated) and the Preparation. The latter refers to the process of creating the necessary files with the partitioning of the mesh depending on the number of processors. These two stages We observe that even though the remeshing and preparation are required just in a few increments they become relatively more expensive than the actual computation because they are performed in sequential. The remeshing takes up to 20% of the cpu time for the periodic meshing and 8% in the automatic case (with 8 processors). The preparation (or partitioning) takes around 40% in both cases.

Figure 4.25 shows at different moments of the simulation how the mesh has evolved due to the displacement of the work piece and the existence of the "Distance from Coil" parameter. For instance, it is observed that after the first 5 increments, most of the zone that was originally covered by the fine mesh is now far from its influence, leaving just a small zone close enough so as to be finely meshed. On the other hand, the automatic remeshing scenario will preserve the fine mesh and distort it. It is just after some 35 increments approximately that the criterion is activated in order to generate the new mesh with less elements. Looking at the evolution of the number of elements in figures 4.26 and 4.27 it is observed that the forced remeshing scenario will decrease the number of elements needed faster than the automatic case (once again depending on the the distance from coil value). This could help reduce the time consumed in the resolution of the Maxwell equations but at the price of introducing the interpolation and transport stages more often.

R-adaptation

One of the advantages of the adaptation algorithm is that, since it is used to artificially move the air mesh and that most of the nodes are dedicated to the discretization of the solids, the loops and iterations are applied only over a reduced amount of nodes. Hence, the cpu time expended on this stage is in any case less than 1% which makes it negligible. In fig. 4.28 and fig. 4.29 is observed that, while the behaviour of the cpu time consumed at each time increment is not the same, it is bounded by an asymptote. This growth is most likely due to the fact that at each increment, the displacement field is initialized to 122 4.8. Computational Time Analysis zero everywhere except for the solid nodes. As the work piece will accelerate due to the growing Lorentz forces, the velocity and thus the displacement per increment will increase. This will translate in more iterations to achieve a satisfactory equilibrium state.

Field Transformation

From fig. 4.23a and fig. 4.24b and table 4.4 we observe that between the interpolation stage and the transformation, it is the former the one that consumes most cpu resource time. The interpolation stage could be largely improved if the remeshing paradigm allowed to preserve the geometrical localization of the partition and refine locally. This way, we wouldn't be forced to research in several (most) of the other partitions for the data.

For the transformation stage we currently transform the magnetic vector potential vector field #» A( #» x ) on the whole domain. Let's remember that the purpose of such operation is to compute the term (σ #» A n ) on the right hand side of Maxwell's equation (previous time step field). One way of optimizing even 123 Chapter 4. Computational Modelling of the Coupled Problem -Objects motion in the immerse FE approach more this stage would be to restrain the operation to only the edges belonging to the work part and the inductor (thus neglecting the air for which σ = 0). Fig. 4.30 -Magnetic vector potential computation at increment 7. Effect of the transport and interpolation after a remeshing step (Forced).

In fig. 4.30 we observe the result of the computation of the magnetic vector potential field #» A at increment 7. A forced remeshing was established from 6 to 7. This means that a mapping was done for field #» A (6) which allowed computing #» A (7) . For the automatic meshing option, the mesh is still the original at this point (no need for mapping). We take as reference the results for the automatic meshing option (auto) computed in sequential (1p). Comparing with the same automatic case but in parallel (2p) a small difference is observed on the field magnitude. The maximum value for the #» A field varies of 0.005%. Which implies that the parallel resolution of Maxwell equations do not pollute excessively the results. On the other hand, we do observe a remarkable difference of -2% for the forced meshing in sequential (1p) and +3% for the parallel forced meshing (2p). This difference is originated on the fact that the meshes are slightly different in their topology after remeshing and also the filtering introduced during the interpolation stage. Such difference, that is cumulated over several increments, will have an impact on the mechanical results as well.

Study Cases

We observe in fig. 4.31 that the displacement field differs on 0.5% between the automatic options in sequential and parallel. For the forced meshing the difference is of 13% for the sequential and of 3% for the parallel option.

Study Cases

Ring Expansion

As in chapter 3 we use the ring expansion case for the validation process. The original geometry presented in Fenton 1996 is used. A squared cross-sectional geometry as shown in fig. 4.32 instead of the circular shaped one used in the previous chapter. Besides, a second case is also studied in which a cylindrical solid wall is included to simulate a collision such as in a joining process.

Free expansion validation

One of the easiest measurable values in an EMF process is the resulting electric current of the system. By means of a Rogowski coil attached to the main circuit it is possible to capture the evolution of this global variable. For this reason we have selected the electric current to be one of the first observables from the simulation. In fig. 4.33a we compare the resulting electric current to the analytic solution. Globally the agreement between both is good, having just a 5% difference on the peak value and 0.5% for the the whole electric current history. Looking at the evolution of the discharge potential at the capacitor level fig. 4.33b, the curves fit almost perfectly. Since the main interest is to characterize the behavior of materials during the EMF process, we now take a look at the thermo-mechanical variables. For this Let us recall that from all the assumptions done for the analytical model, the most important one is the reduction to a 1D problem. All variables are considered uniform throughout the cross section of the solid. This enables focusing only on what happens at the center of the ring. Observing the evolution of the kinematic variables we see that the analytic solution always overestimates the results (see fig. 4.35). This result was already observed in the original work where they developed a 2D axisymmetric code to study the problem. On the other hand, it helps estimating the order of magnitude of the variables and the global behavior.

The advantage of the finite element model resides in its ability to provide us with information on the local phenomena taking place during the forming Chapter 4. Computational Modelling of the Coupled Problem -Objects motion in the immerse FE approach process. When looking at the Lorentz forces for instance, it is found that while there is a major force at the inner section of the ring pushing it forward, there's also a negative force on the opposite side slowing down the deformation, see fig. 4.36. Looking at the stress levels it is observed that at the beginning of the process the stress level is basically homogeneous throughout the cross section (first 10micro seconds). Afterwards, the center and the outer elements undergo higher levels of stress compared with the inner element. fig. 4.36 reveals the lag in the transmission of the body forces. The inner layer (closer to the inductor) will concentrate the magnetic fields due to the skin effect. Nevertheless, the transmission of the wave causes the central and outer region to continue experiencing a higher level of forces while forces in the inner region remain close to zero. The skin depth will also affect the temperature distribution in the body. The Joule heating effect at the skin causes a high increase in temperature locally.

Study Cases

For this case we found 60 • C difference between the outer and the inner nodes.

Another interesting phenomenon revealed in fig. 4.39 is that during the first 5 micro-seconds the Outer temperature is higher than the center one. This is due to the presence of the induced current mainly on the outer surfaces. Then, temperature in the center increases with respect to the outer section due to the transmission of the heat through diffusion, see fig. 4.38.

Constrained expansion: Impact with rigid wall

Considering impact conditions after the rapid expansion allows observing how the deformation load distributes in a more regular fashion along the body.

It was already mentioned in [START_REF] Molinari | Dynamic failure of ductile materials[END_REF] that fast loading helps distributing the necking phenomenon, retarding the damaging and cracking of the material. In fig. 4.40 for the loading case at 7.5kV in free expansion we observe that 2 necks were formed quite near. In quasi-static loading it would have been expected that once a neck has started, the deformation would focus on that single neck. Including the solid wall into the model shows that we reduce the maximun deformation by distributing it more regularly. It will of course depend on several parameters such as the geometry of the work piece and the receptor (wall). Evaluating the radial width compression with respect to the minimal attained value as in fig. 4.41, it is seen that a certain amount of spring back is always present. This spring back has two sources, the elastic recovery of residual stresses and the pressure wave travelling backwards in the impact forming case. For the free forming scenario a decaying oscillating pattern appears, this is a clear sign of an elastic oscillation around an equilibrium point. It is remarkable to see that this equilibrium point rises between 2 and 5kV olts and then drops for 7.5kV olts. It suggests that a transition range exists where the inertia will influence also the small scale elastic oscillations. The impact case evaluated here imposes a rapid stabilization around the final spring back point. This is most likely due to the rapid decay of the mechanical pressure waves.

Another way of analysing the phenomena is by looking at the energy distribution evolution. In fig. 4.42(free expansion column) it is observed that at low input energies (hence, low velocities) the elastic energy is comparable in magnitude with the plastic energy. As the input potential is increased, the elastic energy becomes negligible. This means that most of the energy given to the part has been transformed into irreversible energy-mostly plastic deformation which will also lead to heat dissipation. It is important to see that at the beginning of the process, the input energy brought by Lorentz forces mainly results in an increase of the kinetic energy. Then, the kinetic energy is the responsible of the plastic deformation. This means that the magnetic pressure waves are mainly responsible for the acceleration of the part rather 4.9. Study Cases than directly deforming it. fig. 4.42(impact forming column) shows that at the moment of collision the kinetic energy is suddenly transformed into internal deformation energy, which then mostly ends in plastic deformation.

Magnetic Pulse Welding aluminum on copper (MPW2024)

As was mentioned in the introductory chapter, one of the interests of this technology is the capability for creating high quality joins between dissimilar materials. This can be achieved due to the high speed conditions at impact that change the material behaviour (the solids behave instantaneously like a viscous material). We also have the surface atomic bonding phenomenon which is responsible for the actual joining. Besides, the process happens in such a small span of time that heating effects do not entail phase transformations inside the materials.

The following aims at introducing some of the physical aspects of the process and at the same time some of the difficulties associated with the simulation stage. The model consists in a 5

With Rigid Core

The enormous magnetic fields generated at the skin of the material (82T ⇒ 1.2E6× Earth's field) is strong enough to compress and plastically deform the aluminum tube. The load largely overpasses the yield limit within less than a micro-second (fig. 4.44b). This loading level is yet to be surpassed during impact. According to [START_REF] Zhakhovsky | Two-zone elastic-plastic single shock waves in solids[END_REF] a good estimate for the Hugoniot Elastic limit pressure in aluminum would be P HEL = 12.6Gpa. This value represents the limit between a single phase elastic shock wave and a double phase elastic-plastic wave (see chapter D). It is observed that before impact the material withstands pressures in the range of [0.6 -2.7]Gpa < P HEL . This means that global plastic deformation (before impact) is achieved because the static yield point has been exceeded. Elastic spring back would reappear after the inertia has dissipated and oscillation will happen around the equilibrium point as was observed in the free forming for the ring expansion case. On the onset of impact on the other hand, the pressure reaches levels of > 25GP a (see fig. 4.49). If we assume that P HEL for an Al1050 is similar to the aluminum used in the cited work, then the double shock wave regime has been attained. Meaning that while the elastic shock travels through the material it carries behind a plastic wave that is permanently deforming the body. This plastic shock wave could be the main driving element behind the spring back elimination in high speed forming.

Study Cases

At the end of the process, compressive residual stress in the circumferential direction will remain of about 1GP a contributing to the fixation of the parts. This of course can not be fully evaluated if we don't consider two additional phenomena: (1) the deformation of the solid core and (2) the bonding properties of the surfaces as a function of the contact speed. Changing between a rigid or deformable tool (receptor) has several consequences: From a physical point of view, the impact loading can be large enough as to create severe deformation in the harder body. Consequently, the final shape and quality of the joining section is going to be influenced by this deformation.

From a computational point of view other difficulties show up. First, the fact that it is necessary to compute the deformation field in this new body means an additional cpu requirement. But, most delicate is going to be the management of the contact during impact. In the example shown in fig. 4.45 it can be seen that interpenetration has taken place between the aluminum tube (pink) and the copper tube (green). The cause of this problem is the time step management of the EM module. In the present state of the tool we can only use constant time step in the coupled configuration between the solvers. In FORGE it is necessary to allow automatic time step modification when impact between two deformable bodies is considered. This enables to control the maximum distance of the nodes traveled in a single time step, thus avoiding nodal penetration.

For this example the time step was ∆t = 0.2µs, just before impact the head velocity was about 1.2 × 10 6 mm/s meaning that in the next step we could get a displacement of about 0.25mm. This magnitude is bigger than the mesh size at 4.9. Study Cases the impact zone (0.05mm for the copper and 0.02mm for the aluminum). This excessively large displacement without allowing a modification in the time step is at the center of failure in the computation. One way of bypassing this blocking point is by restarting the thermo-mechanical computation on its own, and using only the velocity computed before impact as the initial condition. This was done as shown in fig. 4.46 where the resulting velocity field as well as the stress, strain, temperature fields were given as initial conditions. The Lorentz body forces and the joule heating were eliminated. From this point forwards the simulation ran only with the thermomechanical module.

] / [ s m V ] stress[GPa principal 1st
Evaluating the Lorentz force at the surface it is observed that the peak level has already bee attained fig. 4.47. Under this condition we can assume that most of the energy has already been transmitted to the work piece. For the proper computation of the mechanics at the impact surface we also made use of a remeshing strategy consisting in defining a "fine mesh size box" fig. 4.48 near the impact zone in order to obtain a small mesh size before impact takes place. As shown in fig. 4.49 this allows to follow the advancing contact front in a detailed manner. In this example a bilateral-sticky contact was defined. This means that once the surfaces have make contact they will remain together. Such approach enables to obtain a first glance at the macro-phenomena and determine the stress levels and deformation suffered by both parts. On the other hand, it overlooks the physics defining the quality of the surface bonding for instance. In [START_REF] Herbst | Suitable Impact Parameters for High-Speed Joining and Influence on the Bonding Zone Microstructure[END_REF] a study is done to correlate the input energy in a EMF process to the quality of the surface welding. The drawback with this correlation is that the input energy can be diffused through several mechanism, so it is difficult to quantitatively know which amount of that energy was actually used to create the macro-bonding. In [START_REF] Wu | Critical velocities for high speed particle deposition in kinetic spraying[END_REF] a study is presented in which the impact velocity between Al-Si particles and a metallic substrate is used as the primary variable to characterize the deposition of the particles. They establish that the creation of the bonding will depend on the competition between the rebound energy and the adhesive energy which are shown in fig. 4.50a. This range also pairs with the gaussian curve for optimal bonding ratio fig. 4.50b.

Returning to the welding case, we saw in fig. 4.46 that the velocities before impact are in the order of [800-1200]m/s. While the results presented by [START_REF] Wu | Critical velocities for high speed particle deposition in kinetic spraying[END_REF] do not represent the conditions of the case at hand, it gives us a direction in the analysis. In their study the optimal bonding speed was found around 700m/s. We can not say that above this point (700m/s) we would encounter the rebounding phenomena by it self in our case. It is important to consider as well the presence of the jet-stream (not seen in the simulation) and the inertia of the matter behind the contact surface. The latter is most probably going to help prevent the rebound adding some pressure whilst impact is in progress. Afterwards, it could play against the bonding by pulling the surfaces apart. Here we point out a possible axis for future studies to generate a joining criteria following a similar scheme as presented before.

Conclusions

theless, since the induction effect drops fast with the gap distance we do not expect a great impact. After the induced force has dropped the main driver of the process is the acquired inertia. The second technique consisted in nodal relocation method based on a non-linear springs analogy. One of the main advantages is that it allows keeping a single mesh through out several time increments with a minimum cost. The time expended solving the relocation of the air nodes is less that 1% of a single increment.

Once the mesh management problem was sorted out, a question that needed to be solved was the transport of the data due to the topological change in the geometrical discretization. The edge data related to the magnetic vector potential posed the issue on how to interpolate values between meshes. The problem was solved by introducing a transport methodology based on three stages: (1) Edge data is transformed into point-wise data on integration points.

(2) the interpolation is done from integration points in the old mesh to integration points in the new mesh. (3) This point data is re-transformed into edge data in the new mesh. This process introduces a new linear system that has to be solved which has the same dimensions as the original electromagnetic problem. Nevertheless, its resolution is less expensive than the original problem since the matrix associated to the transformation is well posed and has a better conditioning.

Another aspect that needs improvement is the interpolation stage. The main issue is linked to the fact that once a new mesh is generated, the parallel partitioning is lost and a new partitioning is created. Since this difference is random, the optimality of data research process is difficult to asses.

In the physical aspect of the process, we saw that the finite element simulation allowed computing the effect of the 3-Dimensional distribution of the magnetic field, which depending on the material and/or the geometry can play against the process. The presence of Lorentz forces on the opposite face of the work part can imply a slow-down on the kinematic acceleration of the part.

It has been stated in literature that the fast loading is at the center of the advantages in forming capabilities of the process. There are two phenomena that play in favour. On one side we have the fast loading itself which transforms into inertial energy. The inertia helps in slowing-down and re-distributing the growth of damage at the microscopic level as has been shown by [Molinari et Chapter 4. Computational Modelling of the Coupled Problem -Objects motion in the immerse FE approach al. 2014]. We gave a closer look at the second mechanism related to the fast constrained arrest of the material. The fast dynamics during impact augments the effect of the inertia, and the constraints (moulds, matrices, dies) help in restraining the deformation of the material.

These phenomena goes hand-to-hand with the reduction of the elastic springback, which is produced by the shock wave mechanics. We could postulate that, if the impact conditions are the appropriate to generate at least a doublewave regime as explained in chapter D, then elastic spring-back elimination can be achieved.

Résumé: Moyennes expérimentales

Ce chapitre est consacré aux installations expérimentales disponibles au laboratoire. La première étape nécessaire afin de modéliser le dispositif expérimental a été sa caractérisation en utilisant le modèle du circuit RLC proposé précédemment. Nous avons vu que cette approche permet une définition simple des paramètres et permet en même temps une identification rapide des paramètres nécessaires. Néanmoins, nous soulignons qu'une différence dans la puissance électrique entre le modèle et le cas réel devrait être prévue pour deux raisons principales: du côté expérimental, si nous ne pouvons capturer qu'un seul pic de courant, il n'y a aucune façon de caractériser le comportement d'amortissement nécessaire à l'identification du paramètre de résistance équivalente. Du point de vue du modèle EF, le fait de tronquer la géométrie de la bobine affecte aussi l'estimation de l'inductance et de la résistance équivalente. Néanmoins, à ce stade, nous tolérons une erreur constatée d'environ 4% sur le courant électrique que nous estimons acceptable pour ce premier modèle.

Les expériences de formage présentées ont été faites dans le but de souligner certains des phénomènes macroscopiques les plus remarquables que l'on observe au cours d'un processus de formage électromagnétique. Le cas de formage indirect a été le modèle livré avec la machine qui permet la formation des matériaux pas ou peu conducteurs (pour l'expérience, nous avons utilisé l'aluminium). Cette configuration présente l'avantage de permettre l'accumulation des énergies d'inerties très importantes dans le poinçon qui formera ensuite la pièce. D'autre part notre intérêt premier était de voir directement l'effet des ondes électromagnétiques sur la pièce de travail. En conséquence, nous avons conçu le cas de gonflage libre en suivant le cas de [Reference]. La simulation du cas de gonflage libre nous fait constater que nous aurons très probablement besoin de réidentifier les données matériau pour les échantillons d'aluminium. Pour le moment, nous sous-estimons la déformation obtenue par rapport au cas expérimental. Cependant, ce résultat pourrait également être amélioré en prenant une section géométrique plus grande de l'inducteur afin de tenir compte du manque de symétrie de révolution de la configuration expérimentale réelle. Enfin, nous avons réalisé un essai d'impression de surface. Il a été observé que l'impact à grande vitesse et l'adoucissement instantané de la matière permet de reproduire des détails -même très petits-de la matrice sur la surface de la pièce.

Chapter 5. Flat sheet forming experiment & tool machining influence

We have discussed how high speed forming of aluminum alloys helps to increase the formability properties of the material as the dynamic nature of the loading changes its mechanical behavior. During impact between the work piece and the die, the material behavior is in a transition between an elasticplastic description and a hydrodynamic one. It means that at the impact zone, the work piece surface can easily adopt not only the global shape of the die but also the small details due to the previous machining of the die. We have witnessed how this phenomenon enables imprinting surface details from a polymeric die into the surface of an A1050 aluminum. We have also tested the imprinting of different surface treatments of a metal matrix into the aluminum sheet.

Experimental setup

EMF Machine

For the experimental settings we have used an electromagnetic forming machine available at the CEMEF laboratory and provided by BMAX industries. The capacitor bank which is currently composed of 12 batteries of 46µF each.

(3) The work station or work bench that is equiped with a pancake-like coil of φ170mm diameter as was introduced in 3.5.

Experimental setup

Identification of the machine's electric parameters

At the beginning of this work we started the development of the simulation software with the aim of piloting the simulation with the same set of parameters that describes (approximately) the electrical behaviour of the real installations. From chapter 3 we saw that the machine can be modelled as shown in fig. 5.2. Identification of C m is easily done from the values reported on the machine. Knowing that all batteries have the same capacity and that they are connected in parallel we can compute the total capacitance by

C m = N i C i = N C = 12 • 46µF = 552µF (5.1)
Identification of R m and L m requires 1. Measure of an electric discharge from the machine in order to determine its actual electric response.

2. Determine a priori how much of the coil is going to be considered in the simulation. As an approximation, R m and L m will contain not only values from the machine but also from the geometrical "truncation" of the coil.

For the measurement of the current output of the machine we perform a discharge without any load, meaning R p and L p = 0. The measurement is carried out with a digital oscilloscope (Handyscope HS4 from TiePie engineering) that is directly connected to a computer for the visualization and processing. Since the magnitude of the current and potential circulating in the main circuit are extremely high a filter is installed in the Rogowski probe1 . Performing a discharge of 2000V olts we obtain the voltage readings shown in fig. 5.3(left). This signal has to be numerically integrated (including the filtering factor) in order to retrieve the actual electric current as shown in fig. 5.3(right).

It has to be remarked that the machine includes a system which automatically disconnects the capacitors from the circuit once they have been fully emptied. This prevents possible damages from backwards currents to the capacitors. Thus the fact that in the potential history plot we see that at 157µs there is a sudden drop in the signal, hence we obtain a single peak of current.

We first estimate the resistance and inductance values of the CAD model of the coil, having previously fixed some values for the machine parameters. From here we get the equivalent values

R eq = R m + R c ; L eq = L m + L c (5.2)
The first guess values for the machine inductance is L m = 69nH. It corresponds to the inductance value reported for the induction bank. Nevertheless,

Experimental setup

there is a large section of conductors not covered in this value, thus the need to recomputed it. For the resistance we do not have any a priory guess value, in that case we simply defined R m = 40mΩ which is close to the one used for the ring expansion case and is also in the same order of magnitude for most reported values in EMF cases.

The perfect harmonic oscillator equation is given by

Q + 2ζω n Q + ω 2 n Q = 0 ; I(t) = Q (5.3)
where

ω n = 1 C m L eq ; ζ = R eq 2 C m L eq (5.4)
and the initial conditions given by

   Q(t = 0) = Q 0 = C m V 0 = 0.1104 Coulombs Q(t = 0) = 0 (5.5)
With this parameters at hand we compute the finite element model and obtain a set of values for the electric current history I F E . The analytic resolution of eq. ( 5.3) produces also a current history data set I analytic . An inverse analysis procedure using the Microsoft Excel optimization toolbox for instance can be done by defining the following objective function

Φ(R c , L c ) = I analytic (R c , L c )dt -I FE dt 2 (5.6)
For which we specify the objective Φ = 0 and the constraints R c > 0, L c > 0.

We then obtain the values From the previous results we validate the FE code not only against an analytical solution from which we find an error of about 0.2%, but also against an experimental result for which the difference is of about 4%. This difference can be attributed to several sources, we point out those we found more relevant.

• The lack of information regarding the decay rate of the oscillations (single peak) makes it difficult to accurately define the actual value of R m .

(Damping factor ζ depends linearly on the resistance).

• The geometrical truncation of the FE model implies a perfect symmetry in the geometry which is not the case for the real model.

• There are several geometrical details that were not included in the CAD design of the coil. Such details have a direct impact on the actual values of R c and L c .

Indirect Forming

As we saw in the introduction of this work, one of the ways of taking advantage of the EMF process is by means of indirect forming. [START_REF] Brower | Electromagnetic devices[END_REF] devised that by placing the part to be formed between a solid metallic punch and an elastomer matrix it could be possible to form several types of materials including non-conductive ones. The first arrangement of the work bench consists on the following mains parts:

• The hammer which is a massive aluminium part responsible of absorbing the electromagnetic energy. Weight: 1087.98gr.

• A steel matrix that has double purpose. A concave zone in the middle for simple forming and four holes in the peripheral for cutting. Weight: 998.10gr.

• A receptor having an external wall of aluminium of 7.5mm thickness and of internal diameter φ74mm. Inside, an elastomer material.

• The work piece that is placed between the matrix and the receptor. For the tests we use φ73mm A1050o sheet of 0.5mm thickness. Weight: 4.6gr It is observed that, for the given discharge conditions, the concave forming at the center of the part is always achieved. The transmitted kinetic energy due to the relatively large mass of the system Hammer+Matrix impacting the part is large enough to guarantee a perfect forming. Additionally, the inertia of the system will cause it to add a compressive load beyond the impact helping to fix the shape. The cutting process on the other hand requires a larger energy in order to be completed. The first reason is due to the material behaviour. In order to cut it is necessary to induce damage through shearing stresses. Full damaging requires, naturally, more energy than simple plastic deformation.

Results

The second reason is of geometrical order. As the diameter of the cutting holes is smaller than the concave section, the resulting force on the section is also smaller.

Free Forming -Bulging experiment

The bulging experiment was inspired from the work of [START_REF] Takatsu | High-speed forming of metal sheets by electromagnetic force[END_REF].

We have carried out several preliminary tests using 160x0.5 mm circular aluminum 1050a sheets for feasibility purposes. After several trials we have found out (results can be observed in fig. 5.10):

• Typical flange wrinkles observed in drawing processes appear for even low discharge potentials.

• The non perfect symmetry of the coil leaves a strong signature on the sheet as can be observed in the irregular formed surface. Anisotropy of

Work Bench

the material due to the previous lamination process may well play an important role here.

• At 2.5kVolts we obtain a peak displacement in the center of the sheet of approximately 3cm. Localized striation is already present in several points around the center.

• For input potentials above 2.6kVolts, we observe the beginning of fracture at the clamping zone due to the section thinning. At 3kVolts full fracture is obtained. Nevertheless, from our preliminary simulation results corresponding to the second bulging scenario we observe in fig. 5.12 that the most critical zones are the border of the clamping zone, where most of the bending load is cumulated and then the center, as also seen in the experimental results. The deformation around the center zone goes through a compression on top/traction on the bottom to a reverse traction on top/compression on the bottom loading. This is caused by the fact that as seen in fig. 5.12a the first zone to move away from the inductor is placed in between the center and the clamping zone (the center is an almost null EM field zone). Thus, at first the center is pulled by the movement of the rest of the part, then, the inertia will force it to keep on advancing. From here the bending regime changes quite fast. This can also be at the origin of the striations that are seen near the center zone in the experimental results. In fig. 5.13 two important results are found. The maximum displacement estimated at the center corresponds to 9.9mm. The experimental case of fig. 5.11 is closer to 14mm. Also, the flange of the sheet did not bend, while in the experiment we expected almost 90 • of bending. The first cause of this difference could be the material data. A new identification would be needed for the used Al sheet. Second, the uneven distribution of the field in the real experiment Chapter 5. Flat sheet forming experiment & tool machining influence that is not taken into account in the simulation.

Work Bench

Surface Imprinting

Note:

The following results are the first findings on a on going research derived from this project and done in collaboration with the University of Northumbria, Newcastle. During the free bulging experiment with the polymeric holder we noticed that the surface patterns of the contact face fig. 5.14 were systematically printed on the metallic work piece whenever an initial gap distance different from zero was left in between, see fig. We have designed a metallic matrix having four kinds of surface finishing in order to quantify the level of details that is possible to print on the work piece.

Work Bench

The different surface treatments are depicted in fig. 5.17. • Turning: Machining process consisting on removal of material through the advance of a rotating tool.

• Milling: Machining process consisting on removal of material through the linear advance of a cutting tool.

• Abrasive blasting: Surface smoothening by blasting a sand jet stream into the surface.

• Grinding: Surface smoothening by friction mechanisms using a rotating tool equipped with a regular distribution of abrasive grains bonded together.

The velocity condition at impact is defined by the kinetic energy given by the Lorentz forces. In our current experimental settings we had no means of measuring directly the velocity or the pressure. The global profile highly depends on the geometry of the coil. For instance, the red path traced in fig. 5.18a shows the zone on the work piece that was above the low-field zones of the coil. It means that this zone had a low kinetic energy and no-impact happened.

We can already observe qualitatively that there is a good printing effect of the turning, milling and abrasive blasting zones on the work flat sheet. On the other hand, the grinding zone cannot be evaluated by the eye because the high quality finishing of the surface makes it almost perfectly flat for the given impact conditions.

Another interesting result is that, as shown in fig. 5.18b for the two lowest quality finishing methods (turning and milling) it is possible to observe their signatures on the rear side of the specimen. This means that the plastic wave was able to travel the 0.5mm of thickness and print the details of the surface finishing.

Conclusions

Conclusions

This chapter was dedicated to the experimental facilities available at the CE-MEF laboratories. The first step needed in order to identify the experimental set-up was the characterization using the RLC circuit model. We have seen that this approach allows a simple definition of the settings and at the same time enables a fast identification of the needed parameters. Nevertheless, we point out that a discrepancy in the electric output between the model and the real case should be expected for two main reasons: From the experimental side, if we only can capture a single current peak there is no way of characterizing the damping behaviour thus this will affect the equivalent resistance parameter. From the FE model point of view, the geometrical truncation of the coil does affect the estimation of the equivalent inductance and resistance. Nevertheless, at this stage we accept the error found of about 4% on the electric current which we estimate acceptable for this first model.

The forming experiments presented so far were done with the purpose of pointing out some of the most remarkable macroscopic phenomena that can be observed during an electromagnetic forming process. The indirect forming case was the model delivered with the machine that enables forming nonconductive materials (for the experiment we used aluminium). This set-up has the advantage of allowing accumulating large inertial energy in the punch part that will form the work piece. On the other hand our first interest was to see the effect of the electromagnetic waves on the work piece directly. Hence we designed the free bulging case following [START_REF] Takatsu | High-speed forming of metal sheets by electromagnetic force[END_REF]. From the simulation of the free bulging case we saw that we most likely need to reevaluate the material data for the aluminium samples. For the moment we underestimate the deformation obtained with respect to the experimental case. However, this result could also be improved by taking a larger geometrical section of the inductor in order to take into account the lack of axisymmetry of the actual experimental configuration. Finally, we found quite interesting the results of the surface printing experiment. It was observed that the high-speed impact and the instantaneous material softening allows copying even small details from the die/matrix into the work piece.

Chapter 6

General conclusions and perspectives

This Ph.D. work had as main objective to develop a predictive computational tool capable of simulating the electromagnetic forming process taking into account all the complexities of multiphysical interactions within a 3-Dimensional framework. This was motivated from the fact that while we know from literature and experience that the process is promising for the metal forming industry, its multiphysical nature makes it complex to design and understand why it works so well.

Two key points are at the center of the discussion:

(1) The increased forming limits. This phenomenon depends on the competition between inertial forces and damage propagation. Since the mechanics of damage have to be studied at the micro-scale we left this subject out of our current model. Nevertheless, we can point out as perspective to Perform a multi-scale study of the process in which the damage phenomenon is modelled with the representative volume element approach (RVE) and link this information to the macro-scale model. A team in our laboratory is currently working on the development of the tool Digi-µ for meso-scale modelling of materials which could be linked to our softwares and could serve to the purpose.

(2) The spring-back elimination: This phenomenon is linked to the high-speed impact of the solid material with the die. The first driving factor is the propagation of the plastic shock-wave and the second one is the creation of compressive stresses helping in retaining the parts together (we should add also Chapter 6. General conclusions and perspectives the surface bonding for cases where bonding conditions are achieved).

Perspective: To perform a fast dynamic simulation at the nano-second scale to accurately witness the effect of Lorentz forces propagating from one side of the material, the impact shock wave coming from the rear side and both contributing to the formation of the residual compressive stresses. This study could be done with the tool developed in this work by carrying out the simulation with the appropriate time scale.

Since the main focus of the project was the numerical development of the simulation tool, most of the attention was given to the mathematical and computational aspects. In chapter 2 we reviewed the fundamental principles behind the description of the electromagnetism and the thermo-mechanics. We left a question open regarding the derivation of Maxwell's equations: Are they the actual equations that should be used in a continuum description of deformable objects? We proposed a continuum derivation of the electromagnetic equations following what is known as the Heaviside-Hertz model. From it we were able to derive the Lorentz force on a particle as a result of the fields equations without additional terms. In [Christov 2006] was mentioned that the Lorentz contraction could be explained starting from the Heaviside-Hertz model. In the appendix B.2 we show an alternative approach to derive the contraction factor using this model.

Besides the fundamental discussion rising from the use of the model we also found that in the domain of magneto-hydrodynamics the electric potential equations is usually presented as ∇ • (∇φ) = ∇ • ( #» v × #» B) which can also be obtained from the Heaviside-Hertz model as done in appendix B.1. This raises a new question (perspective): Are there EMF configurations in which the velocity of the deformable object generates an induced electric potential? If so, we should solve eq. (B.2) on the inductor and the work piece as well.

In chapter 3 we established the finite elements framework for the simulation of the eletromagnetic problem using a fully immersed approach. We also dedicated a section to the modelling of the electric machinery behind the process with the aim of piloting the simulation in the same way we pilot the real experiment. The resolution of the coupled Finite Elements-Circuit equation is done by means of a fixed-point scheme. We see as a perspective for generalization of our code to modify the resolution by a constrained resolution of the linear system, which would enable to solve simultaneously the vector potential ( #» A) and the electric current in the inductor (I).

From the comparison with the electric results of the experimental case we saw that defining the correct set of parameters for the circuit equation is not trivial because we need to take into account the geometrical simplifications done in the CAD model with respect to the actual coil. Nonetheless, the identification can be done by using an analytical harmonic oscillator model, from where we can get a good initial approximation of the values.

Coupling between the electromagnetic problem and the thermo-mechanical problem was introduced in chapter 4. We selected a sequential coupling strategy based on independent mesh management for the different physics. This created the need for a mesh management method allowing the mesh of the air, in the EM computation, to track the displacement of the solids. Two strategies were implemented. An r-adaptation technique enabling to define a mesh nodal displacement and a remeshing that generates a new mesh when the virtual displacement creates excessive distortion.

Perspective:

The remeshing algorithm should be rewritten in order to optimally work in parallel. It currently constitutes a bottleneck for the computation because it is done sequentially. Then,  

 if i = l ⇒ ∂x i ∂x l = 1 if i = l ⇒ ∂x i ∂x l = 0 (A.6)
For any case ∂ ∂x 1 ∂x i ∂x l = 0 meaning that ∂ ∂x l J êl = 0, from where it follows that

d dt J = ∂ ∂t J (A.7)
For the next part of the demonstration we will make use of the Jacobi formula In the domain of magneto-hydrodynamics we can find quite often the use of the following expression for the electric potential (see [Moreau 1990;[START_REF] Shatrov | Three-dimensional linear stability analysis of lid-driven magnetohydrodynamic cavity flow[END_REF][START_REF] Jin | Three-dimensional Flow in a Driven Cavity Subjected to an External Magnetic Field[END_REF] for instance)

∇ • (∇φ) = ∇ • ( #» v × #» B) (B.1)
From the proposed development in chapter 2 we can easily obtain it by taking the divergence of the Maxwell-Ampere's Law

∇ • ∇ × #» H = ∇ • ∂ #» D ∂t -∇ × ( #» v × #» D) + #» J 0 = ∂ ∂t ρ e + ∇ • ( #» J ) 0 = ∂ ∂t ρ e -∇ • (σ∇φ) -∇ • (σ∂ t #» A) + ∇ • (σ #» v × ∇ × #» A)
From where ∂ ∂t ρ e is usually set to zero and, by means of the coulomb gauge we expect ∇ • (σ∂ t #» A) = 0 as well, thus finding We obtain

∇ • (σ∇φ) = ∇ • (σ #» v × ∇ × #» A) (B.
1 + (ω + k • v) 2 c 2 -k2 -1 ( k - 1 c 2 (2ω + k • v)v) • k = 0  k 2 - ω 2 + k • v 2 + 2ω k • v c 2   -1 k - 1 c 2 2ω k • v + k • v 2 = 1  k 2 - ω 2 + k • v 2 c 2   -1 k - 1 c 2 k • v 2 = 1 1 k2 -1 c 2 ω 2 + ( k • v) 2 kT Ī - 1 c 2 v ⊗ v k = 1 (B.3)
Imagine the particular case of v = v, 0, 0 T , then

1 = 1 k2 -1 c 2 ω 2 + ( k • v) 2    1 - v c 2 k x 2 + k 2 y + k 2 z    (B.4)
Equation (B.3) describes the geometry of the propagation wavefront around and observer travelling with the source as described from the reference or 'stationary' frame. We can obtain from it the angles of the main propagation axes.

In the particular case of (B.4) it means that the wavefront propagation is an ellipsoid contracted in the direction of movement of the source by the factor

1 -v c
2 . Given Ē = E 0 n 1 , n 2 , n 3 T ; n i = cos(φ i ) ; i = 1 . . . 3

Then

n 1 = 1 -v c 2 k x k2 -1 c 2 ω 2 + ( k • v) 2 n 2 = k y k2 -1 c 2 ω 2 + ( k • v) 2 n 3 = k z k2 -1 c 2 ω 2 + ( k • v) 2 (B.5)

D.1 Euler equations

In chapter 2, when developing the governing equations of the thermo-mechanical problem, at all times we presented the equations in its Lagrangian and Eulerian forms. The purpose was to set up the background framework from which we can then develop the rest of the phenomena. In this case we are interested in the shock wave propagation due to impact conditions. The study of propagating waves is done easily when regarding them as moving entities within a fixed frame as in fluid dynamics. For this reason it is convenient to use the Eulerian representation of the governing equations. The key concept was the introduction of the jump condition in 1D for a conserved quantity D Dt Where ρ represents the density of the material in the shock affected zone (1) or the non affected zone (2). u is the velocity of a particle in each of the zones.

Assuming that the non-affected zone is in a static state with respect to the affected zone we could say u 2 = 0, ρ 2 = ρ 0 , ρ 1 = ρ and call u 1 = u p (p is the standard notation for the particle velocity). It follows that the compression can be measured by ρ 0 ρ = u su p u s (D.13)

The Hugoniot is an important relation in the construction of the EOS. At moderate pressures it can be represented by the linear relationship [START_REF] Kanel | Shock-wave phenomena and the properties of condensed matter[END_REF].

u s = c + su p (D.14)
Where c is some1 measure of the speed of sound and s is a proportionality constant which has a value between 1 and 1.7 for most of the known materials. Another measure of the compression is given by the variable

χ = 1 - ρ 0 ρ = u p c + su p = 1 s + c up (D.15)

D.3 The Pv diagram for shock mechanics

Determining the state diagram of the shock process enables to understand the different phases through which the material goes when submitted to shock loading. In the works of [START_REF] Zhakhovsky | Elastic-plastic phenomena in ultrashort shock waves[END_REF][START_REF] Zhakhovsky | Two-zone elastic-plastic single shock waves in solids[END_REF] [ [START_REF] Zhakhovsky | Elastic-plastic phenomena in ultrashort shock waves[END_REF][START_REF] Zhakhovsky | Two-zone elastic-plastic single shock waves in solids[END_REF] HEL: Hugoniot Elastic Limit EZ:

Elastically compressed metastable elastic state HEL * : Hugoniot final metastable elastic limit S el : Hugoniot Elastic Adiabat S pl : Hugoniot Plastic Adiabat OD: Over driven OD * : Over driven limit of the two-phase shock wave D.3. The Pv diagram for shock mechanics This diagram resumes a large amount of work dedicated to the identification of the shock regimes that starts from a single elastic wave, then goes to a elastic-plastic splited regime, and then a single plastic wave. The last is still under debate. Shock Regimes • Gradient ∇ϕ i = F∇ η,γ,ξ ϕ i (E.4)

• Vector interpolation functions

Ψ 1 = ϕ 1 ∇ϕ 2 -ϕ 2 ∇ϕ 1 Ψ 2 = ϕ 2 ∇ϕ 3 -ϕ 3 ∇ϕ 2 Ψ 3 = ϕ 3 ∇ϕ 1 -ϕ 1 ∇ϕ 3 Ψ 4 = ϕ 1 ∇ϕ 4 -ϕ 4 ∇ϕ 1 Ψ 5 = ϕ 2 ∇ϕ 4 -ϕ 4 ∇ϕ 2 Ψ 6 = ϕ 3 ∇ϕ 4 -ϕ 4 ∇ϕ 3 (E.5) • Curl ∇ × Ψ d = -2∇ϕ i × ∇ϕ j (E.6)
Then, the local element matrices defined in 3.33 are

M 11 = 1 30α 2 h 2           1 -1 4 -1 4 0 0 0 α 2 +1 2 -α 2 4 -α 2 -1 4 -α 2 -1 2 -α 2 -1 2 α 2 +1 2 α 2 -1 2 α 2 -1 4 α 2 -1 2 α 2 +4 2 α 2 +4 4 2α 2 +3 4 α 2 +4 2 2α 2 +3 4 3α 2 +2 2           (E.7) eig(M 11 ) = 1 30α 2 h 2           α 2 2 + 3
M cc = 2 3α 4 h 4           1 0 0 -1 1 0 α 2 0 0 -α 2 α 2 α 2 α 2 0 -α 2 α 2 + 1 -1 -α 2 α 2 + 1 -α 2 2α 2           (E.9) eig(M cc ) = 2 3α 4 h 4           α 2 - √ 4α 4 -4α 2 +9 2 + 3 2 α 2 + √ 4α 4 -4α 2 +9 2 + 3 2 0 0 0 4α 2           (E.10) M 10 = 1 24α 2 h 2           -1 1 0 0 0 -1 α 2 -(α 2 -1) 1 0 -1 α 2 -1 -2 -1 -1 α 2 + 3 -1 -2 -1 α 2 + 3 -1 -1 -2 2(α 2 + 1)           (E.11)

E.1 Characteristic Matrices

Note: For the sake of clarity, only the non-zero elements of the matrices are listed below.

Fig. 1

 1 Fig. 1.1 -Some of the first applications to the electromagnetic forming applications according to Zittel 2010.

Fig. 1 . 2 -

 12 Fig. 1.2 -Fundamental principle behind high speed forming processes.

Fig. 1 . 3 -

 13 Fig. 1.3 -Under proper impact velocity conditions a jet stream appear that cleans the contact zone.

Fig. 1 . 4 -

 14 Fig. 1.4 -Classic forming configurations as described by Harvey et al. 1961.

Fig. 1 . 5 -

 15 Fig. 1.5 -Stages of the EMF forming process.

Fig. 1 . 6 -

 16 Fig. 1.6 -Example forming diagram for a 6061 T4 aluminum as taken fromBalanethiram et al. 1994.

Fig. 2 . 1 -

 21 Fig. 2.1 -Decomposition of the main elements in the electromagnetic forming process.(Left) The EMF machine (Command+energy storage). (Center) Work bench (Coil+matrix). (Right) Work piece after deformation and strain map.

Chapter 2 .

 2 Physical & Mathematical Frameworkthe geometrical evolution since the intensity of the EM fields decreases rapidly with distance. This interaction is summarized in fig.2.2.

Fig. 2

 2 Fig. 2.2 -Interaction of the main different physics involved in the EMF process

Fig. 2 . 3 -

 23 Fig. 2.3 -Typical reference systems in continuum mechanical description. The Eulerian or static frame of reference and the Lagrangian or deforming frame.

  Magnetic flux intensity ρ e : Electric charge density Note that we have eliminated the term associated to the current density in (2.15)b. This term will re-appear when we properly develop the time derivative of the integral.

4 -

 4 Fig. 2.4 -Additive strain model.

σ a :

 a Theoretical flow stress at 0 Kelvins K = 0 for fcc metals n : Strain hardening exponent µ T , µ 0 : Shear modulus at temperatures T and reference σ t : Path-dependent part φ : Empirical state-dependent part Σ : Previous loading stress state

Fig. 2

 2 Fig. 2.5 -Interaction of the main different physics involved in the EMF process. Mathematical representation.

Fig. 3

 3 Fig. 3.1 -Boundaries of a EMF process.Ω represents the global domain solids+surroundings. Ω p is the work piece. Ω I represents the inductor domain. The electrical input and output connections of the inductor are given by Γ I inp and Γ I out . After[START_REF] Biro | On the use of the magnetic vector potential in the finite-element analysis of three-dimensional eddy currents[END_REF]].

Fig. 3

 3 Fig. 3.3 -2D triangular element.

Fig. 3 . 4 -

 34 Fig. 3.4 -Edge base functions in a 2D P1 triangular element.

Fig. 3 . 5 -

 35 Fig. 3.5 -Electrical scheme of a general EMF process

  Fig. 3.6 -Ring expansion model configuration.

Fig. 3

 3 Fig. 3.7 -Finite element model. 15 • section approximation.

Fig. 3

 3 Fig. 3.8 -Different refinement levels proposed. Global mesh size of 4mm. Local mesh size: 0.1mm (top, M1), 0.05mm (center, M2), 0.025mm (bottom, M3).

  L m = 115nH.

Fig. 3

 3 Fig.3.9 -Number of iterations in the first increment. Influence of χ (for fixed θ = 0.5) and θ (for fixed χ = 0.5).

  Fig. 3.10 -Evolution of the electric current. Influence of the Time Step. Machine inductance L = 0nH.

Fig. 3 .Fig. 3 .

 33 Fig. 3.12 -Evolution of the electric potential. Influence of the Time Step. Machine inductance L = 0nH.

  Fig. 3.14 -Analytic energy balance at the inductor-work piece level. Case1.

  Fig. 3.15 -Analytic energy balance at the inductor-work piece level. Case2.

  Fig. 3.16 -Energy balance from numerical time integration of the current and potential outputs. Case 1. Straight lines correspond to the analytical solution, Mark plots correspond to the FE output.

  Fig.3.17 -Energy balance from numerical time integration of the current and potential outputs. Case 2. Straight lines correspond to the analytical solution, Mark plots correspond to the FE output.

Fig. 3 .

 3 Fig. 3.18 -Scheme of the domain partitioning in the parallel computation context. The meshing program generates sub-domains in which elements belong to a single domain, nodes are shared at the interface and edges are also shared following nodal partitioning paradigm.

Fig. 3 .

 3 Fig. 3.19 -Speed up of parallel computing.

Fig. 3 .

 3 Fig. 3.21 -Model mesh.

Chapter 4 .

 4 Computational Modelling of the Coupled Problem -Objects motion in the immerse FE approach forces such as gravity and Lorentz induced forces.

Fig. 4 .

 4 Fig. 4.1 -Boundary conditions of the mechanical problem

Fig. 4 .

 4 Fig. 4.2 -Scheme of the penetration distance condition between colliding objects.

Fig. 4 .

 4 Fig. 4.3 -Finite elements discretization of a physical domain Ω into smaller subdomains Ω k .

Fig. 4 . 4 -

 44 Fig. 4.4 -2D scheme of a P +1 /P 1 element. The velocity is defined for all the nodes plus the central node (bubble) within the element (P + 1 ) and the pressure is defined only for the nodes (P 1 ).

Fig. 4 .

 4 Fig. 4.5 -Iterative time evolution scheme for implicit resolutions.

Fig. 4 .

 4 Fig. 4.6 -Boundary conditions of the heat transfer problem.

Figure 4 .

 4 Figure 4.8 shows the global work flow to accomplish a single time step increment. The coupling between the solvers is driven by the communication of: The Lorentz volumetric forces ( #» F L ) and the Joule heating term ( QJ ) at each element of the solid body from the EM module. The temperature evolution (T ) and the velocity field ( #» v ) or the new position ( #» x ) of the nodes of the solid body are transmitted from the TM module.

Fig. 4 .

 4 Fig. 4.8 -Global algorithm of the sequential coupling algorithm.

Figure 4 .

 4 Figure 4.9 shows an example schematizing the actual configuration of the finite element modeling. The electromagnetic module handles the electric source (inductor), the work piece and the air for the transmission of the EM waves. In the thermo-mechanical module we only need to define the work piece if it corresponds to a free forming case and if necessary any additional object that serves as tool or matrix for the final shape.

Fig. 4 .

 4 Fig. 4.9 -Configuration of the finite element model for the MPW2024 joining case.

Chapter 4 .

 4 Computational Modelling of the Coupled Problem -Objects motion in the immerse FE approach

Fig. 4 .

 4 Fig. 4.10 -Skin depth concentration of the electric current.

Fig. 4 .

 4 Fig. 4.11 -Automatic mesh generation by material parameters considerations. In blue the air mesh, in red the coil and green the work piece.

Fig. 4 .

 4 Fig. 4.12 -Stages of the nodal relocation method. The surrounding air mesh Ω has to follow the movement of the internal solids.

Fig. 4 .

 4 Fig. 4.13 -Springs based method.

Fig. 4 .

 4 Fig. 4.14 -Scheme of boundary node between two neighboring partitions domains of the mesh.

Fig. 4 .

 4 Fig. 4.15 -Symmetry plane constraint.

Fig. 4 .

 4 Fig. 4.16 -A given field has to be transferred between two coincident but nonidentical meshes.

  .

  (a) Field Transformation: from edge to point data (b) Field transport or recovery (mesh to mesh), and (c) Field Transformation: from point data to edge data.

Chapter 4 .

 4 Fig.4.17 -Stages of the mesh-to-mesh data transfer process.

( a )

 a Reference element in 3D (b) Scheme of the overlap.

Fig. 4 .

 4 Fig. 4.18 -Reference element for the interpolation of a given field from the integration points.

Fig. 4 .

 4 Fig. 4.19 -2D 2 elements test for the field transformation methods.

Chapter 4 .Fig. 4 .

 44 Fig. 4.20 -Modification of the domain partitioning among the different processes between two time-step increments. On the left a first hypothetical partitioning with its circumscribing circles (spheres in 3D) at increment n and to the right a new partitioning at n + 1.

Fig. 4 .

 4 Fig. 4.21 -Estimation of the overlapping between two circumscribing spheres.

Fig. 4 .

 4 Fig. 4.22 -Partitioning of a 2D mesh in four domains and circumscribing spheres. (left) First partitioning, (right) Second partitioning.

  75 0.25 0.35 0.33 0.50 0.60 4.8. Computational Time AnalysisTable 4.2 -Comparison of the second partitioning with the original by means of distance (D), coverage ratio (CR) and the resemblance factor (RF).

  (a) Forced remeshing every 10 increments and (b) automatic remeshing by verification of a quality criterion defined by the rule eq. (4.64). Each case is compared in sequential and in parallel. The main comparison variable is going to be the total (fig. 4.23b,fig. 4.24b) and average (fig. 4.23a,fig. 4.24a) CPU time expended in each stage.

Fig. 4 .

 4 Fig. 4.23 -Time consumed by each stage of the simulation for any increment in the case of forced remeshing.

Fig. 4 .

 4 Fig. 4.24 -Time consumed by each stage of the simulation for any increment in the case of automatic remeshing.

Chapter 4 .Fig. 4 .

 44 Fig. 4.25 -Mesh adaptation to the solid body displacement during forced remeshing (top) and automatic remeshing (bottom).

Fig. 4 .Fig. 4 .

 44 Fig. 4.26 -Evolution of the number of elements in the forced remeshing scenario.

Fig. 4 .Fig. 4 .

 44 Fig. 4.28 -Evolution of the CPU time consumed for the r-adaptation algorithm in the forced remeshing scenario.

Fig. 4 .

 4 Fig. 4.31 -Displacement field of the work piece at the end of increment 100. Comparison of the effect of the electromagnetic field computed under automatic or forced remeshing conditions.

Fig. 4 .

 4 Fig. 4.32 -Configuration of the finite element coupled models for the ring expansion case.

Fig. 4 .

 4 Fig. 4.33 -Comparison between the analytical model and the 3D FEM. Electrical

Fig. 4 .

 4 Fig. 4.34 -Selection of output positions. Inner, Central and Outer labels used for identification of the the variables with respect the ring radius.

Fig. 4 .

 4 Fig. 4.35 -Comparison between the analytical model and the 3D FEM. Kinematics.

Fig. 4 .Fig. 4 .

 44 Fig. 4.36 -Evolution of the Lorentz body forces. Comparison Analytic vs. FEM.

Fig. 4 .

 4 Fig. 4.38 -Evolution of the Joule heating term. Comparison Analytic vs. FEM.

Fig. 4 .

 4 Fig. 4.40 -Comparison of the equivalent plastic deformation in constrained (top) vs. free (bottom) forming.

Chapter 4 .Fig. 4 .

 44 Fig. 4.41 -Width variation after minimal compression. (Left) Free expansion (Right) Impact forming.

Fig. 4 .Chapter 4 .

 44 Fig. 4.42 -Energy density balance evolution. From top to bottom, at 2kV, 5kV, 7.5kV and 10kV

Fig. 4 .

 4 Fig. 4.43 -Joint of a aluminum A2024 with a copper tube due to high speed loading.

( a )Fig. 4 .

 a4 Fig. 4.44 -Magnetic and mechanical loads on the aluminum tube.

Fig. 4 .

 4 Fig. 4.45 -Mesh distortion during impact between the two deformable bodies when no automatic time-stepping is allowed. EM model (top) TM model (bottom).

Fig. 4 .

 4 Fig. 4.46 -Thermomechanical simulation re-started in decoupled fashion just before impact to allow automatic time-stepping.

Fig. 4 .

 4 Fig. 4.47 -Evolution of the Lorentz force density as computed in a surface element of the aluminum tube. Impact at t = 2.4µs (red line).

Fig. 4 .

 4 Fig.4.48 -Fine mesh box (pink) for the refinement of the close-to-impact zone in the aluminum tube.

Fig. 4 .

 4 Fig. 4.49 -Evolution of the contact zone during impact.

  Machine . . . . . . . . . . . . . . . . . . . . . . . . 5.1.2 Identification of the machine's electric parameters . . . 5.2 Work Bench . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Indirect Forming . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Free Forming -Bulging experiment . . . . . . . . . . . . 5.2.3 Surface Imprinting . . . . . . . . . . . . . . . . . . . . . 5.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Fig. 5 . 1 -

 51 Fig. 5.1 -Electromagnetic forming machine

Fig. 5

 5 Fig. 5.2 -Electrical scheme of a general EMF process

  Fig. 5.3 -Electric Output from the EMF machine in load free test at 2000V olts.

  Fig. 5.6 -Comparison of electric current experimental vs. analytical and FEM.

  Fig. 5.8 -Indirect forming work bench.

Fig. 5

 5 Fig. 5.9 -Indirect forming test results.

2014CEMEF

  Fig. 5.10 -Experimental results of a flat sheep metal expansion without free flange. Trails at 1.5kVolts (top-left), 2.5kVolts (top-right), 2.65kVolts (bottom-left) and 3kVolts(bottom right).

Fig. 5 .

 5 Fig.5.12 -Deformation under the free bulging case. Load 2.5kV olts

Fig. 5 .

 5 Fig. 5.13 -Maximum displacement. Load 2.5kV olts. (Expanded view of the 2 • model)

Fig. 5 .

 5 Fig. 5.14 -Polymer holder for free expansion testing (fig. 5.14a). Surface profile of the contact zone (fig. 5.14b).

  Fig. 5.15 -Work piece after free bulging fig. 5.15a and measurement of visible striations on the surface of the work piece fig. 5.15b.

Fig. 5 .

 5 Fig. 5.16 -Surface roughness profile of the polymeric holder.

Fig. 5 .

 5 Fig. 5.17 -Surface imprinting matrix. Description of the different surface finishing.

Fig. 5 .

 5 Fig. 5.18 -Deformed sheet at 3.5kV and with a 2mm gap.

  t) is a given tensor field dependent on t and Adj A is the adjoint matrix of A. Then,

  + ∇(ρv 2 ) = ∇ • S -∇p + ρ b (D.2) Energy conservation ∂ ∂t (ρe) + ∇ • [v (ρe + p)] = ∇ • (S • v) + ρv • b + ρ Q -∇ • q (D.3)The well known Euler equations for fluid dynamics are obtained under the following hypothesis:When the theory was first established, all the mathematical results were obtained mainly in 1D approximations. This allowed to generate analytic equa-D.2. Rankine-Hugoniot shock theory tions and properly formulate the fundamental concepts. The first step was to simplify the Euler equations into a 1D form:

  Fig. D.1 -Scheme of a two-zone shock wave propagating within a solid material.

  we have found one of the most complete descriptions as seen in fig. D.2.

  Fig. D.2 -Schematic of Pv shock Hugoniot with elastic and plastic fronts.[[START_REF] Zhakhovsky | Elastic-plastic phenomena in ultrashort shock waves[END_REF][START_REF] Zhakhovsky | Two-zone elastic-plastic single shock waves in solids[END_REF]

  Figure fig. D.3 is an example state diagram for an aluminium sample comparing experimental and numerical data.

  Fig. D.3 -(From[START_REF] Zhakhovsky | Two-zone elastic-plastic single shock waves in solids[END_REF]): P-V Hugoniot for perfect AL[111] crystal and a sample with a vacancy concentration of 10 -3 .

4 B 4 M 4 M

 444 22 = 1 B 23 = -b 1 B 33 = ∆tµ 0 ε 0 ε r m 11 b 4 + σm 11 b 1 + m cc µ r B 35 = ∆tµ 0 m 10 G σ + ε 21 = b 6 M 22 = b 5 M 23 = -b 31 = -m 11 ∆tµ 0 (b 3 σ + b 6 ε 0 ε r ) M 32 = -m 11 ∆tµ 0 (b 2 σ + b 5 ε 0 ε r ) M 33 = m 11 ∆tµ 0 (b 1 σ + b 4 ε 0 ε r ) M 35 = ∆tµ 0 m 10 G

4 The Solid Mechanics and Heat Transfer Problems . . . . . .

	Chapter 2
	Physical & Mathematical
	Framework

sense, or not converge at all. If they do, we know that the solution does correspond to the equilibrium state Contents 2.1 A coupled multiphysics problem . . . . . . . . . . . . . . . . 2.2 Continuum Mathematics . . . . . . . . . . . . . . . . . . . . . 2.2.1 Reference frames . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Reynold's Transport theorem -Volume Integral time derivative . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 The Electromagnetic Problem . . . . . . . . . . . . . . . . . . 2.3.1 Electromagnetic Equations . . . . . . . . . . . . . . . . 2.3.2 Classical electromagnetism . . . . . . . . . . . . . . . . 2.2.4.1 Mass conservation: . . . . . . . . . . . . . . . . . . . . . 2.4.2 Linear momentum conservation: . . . . . . . . . . . . . 2.4.3 Angular momentum conservation: . . . . . . . . . . . . 2.4.4 Energy conservation, 1st Thermodynamic Principle: . . 2.4.5 Material physical modelling . . . . . . . . . . . . . . . . 2.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 2. Physical & Mathematical Framework

. The Solid Mechanics and Heat Transfer Problems Definition 1. Bulk

	2.4Modulus (κ): Represents the resistance to uniform compres-
	sion	κ = ρ	dp dρ	(2.61)
	Where p represents the pressure felt by the differential volume element. This
	definition is used to track the evolution of the volume by insertion of eq. (2.61)
	into eq. (2.35) from where
		tr σ I	(2.60)
	Compressibility In solids, under elastic deformation it is possible to uni-
	formly compress a volume element without modifying its shape. This prop-
	erty is measured by

. The Solid Mechanics and Heat Transfer Problems energy

  and suppose that it depends linearly on the temperature

	2.4
	• Energy evolution d dt u int : We neglect any chemical transformation of the

Chapter 3 Computational Modelling of Electromagnetism Contents 3.1 Electromagnetic Modelling . . . . . . . . . . . . . . . . . . .

	3.1.1 The potential formulation . . . . . . . . . . . . . . . . .
	3.1.2 Boundary conditions . . . . . . . . . . . . . . . . . . . .
	3.1.3 Finite elements for electromagnetism . . . . . . . . . .
	3.1.4 Time integration schemes . . . . . . . . . . . . . . . . .
	3.2

Machinery modelling and implementation . . . . . . . . . . 3.3 Validation of the EM finite elements module . . . . . . . . .

	3.3.1 Description . . . . . . . . . . . . . . . . . . . . . . . . .
	3.3.2 1D Analytical Model . . . . . . . . . . . . . . . . . . . .
	3.3.3 3D Finite Elements Model . . . . . . . . . . . . . . . . .
	3.

4 Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . .

	3.4.1 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	3.4.2 Energy Analysis . . . . . . . . . . . . . . . . . . . . . . .
	3.5

Parallel computation . . . . . . . . . . . . . . . . . . . . . . .

	3.5.1 Partitioning paradigm . . . . . . . . . . . . . . . . . . .
	3.5.2 Test case . . . . . . . . . . . . . . . . . . . . . . . . . . .
	3.6

Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 3. Computational Modelling of Electromagnetism Résumé: Modélisation numérique du phénomène élec- tromagnétique

	Dans ce chapitre, nous examinons et mettons en place les outils numériques
	nécessaires pour modéliser le problème électromagnétique par la méthode
	des éléments finis. En raison de la présence de l'opérateur rotationnel dans les
	équations de l'électromagnétisme, il est nécessaire d'adapter la formulation
	des éléments finis, plus précisément à l'aide des éléments d'arête de Nedelec
	Nédélec 1986. Ces éléments éliminent les solutions parasites qui surgiraient
	si des éléments finis nodaux classiques étaient utilisés Jin 2002 et simplifient
	également la définition des conditions aux limites électromagnétiques.

Table 3 .1 -Process parameters Machine

 3

Table 3 .

 3

	2 -Ring expansion case geometric
	parameters
	Geometry
	Coil mean radius a 1 [mm]	14
	Coil cross-section radius b 1 [mm] 0.5
	Ring mean radius a 2 [mm]	16
	Ring cross-section radius b 2 [mm] 0.5

Table 3 .

 3

		Nodes Elements	Edges
	M1	7535	42109	50390
	M2 101564	598278	701895
	M3 457931 2711134 3175352

3 -Stats of the FEM models.

Table 3 .

 3 4 -Current Peak deviation from the analytic model. Time step influence.

Table 3 .

 3

	5 -Eigenvalues of the integration	Table 3.6 -Eigenvalues of the integration
	methods. L = 0nH			methods. L = 115nH
		QSEI QSSC DYN		QSEI QSSC DYN
		-	-	-1		-	-	-1
	FEM	0.164 0.303 -0.436	FEM	0.164 0.303 -0.436
		0	0	0		0	0	0
	Circuit	-0.034 -0.034 -0.034 0 0 0	Circuit	0.923 0.923 0.923 0 0 0

Chapter 4. Computational Modelling of the Coupled Problem -Objects motion in the immerse FE approach géométrique

  . Les données liées au potentiel-vecteur magnétique posent la question de l'interpolation des valeurs entre les maillages. Le problème a été résolu par l'introduction d'une méthode de transport basé sur trois étapes:1. Les données aux arrêtes se transforment en données de points sur les points d'intégration.2. l'interpolation se fait entre les points d'intégration de l'ancien maillage et les points d'intégration du nouveau maillage.3. Ces données aux points d'intégrations sont transformées en données aux arêtes dans le nouveau maillage.Ce procédé introduit un nouveau système linéaire qui doit être résolu et qui a les mêmes dimensions que le problème électromagnétique d'origine. Néanmoins, sa résolution est moins chère que le problème initial puisque la matrice associée à la transformation a un meilleur conditionnement.

Un autre aspect qui a besoin d'être amélioré est l'étape d'interpolation. Le problème principal est lié au fait qu'une fois qu'un nouveau maillage est généré, le partitionnement parallèle est perdu et une nouvelle partition est créée. Étant donné que ce partitionnement est aléatoire, l'optimalité du processus de recherche des données est difficile à évaluer.

Table 4

 4

	.3 -Total cpu times of the electromagnetic computation stages in coupled sim-
	ulation.
	Processors	1	2	4	8
	Forced	21'24" 12'26" 09'37" 10'11"
	Automatic 17'32" 10'05" 06'49" 05'34"

4.8. Computational Time Analysis are

  currently carried out sequentially, thus they are almost independent on the number of partitions. They are mainly affected by the size of the resulting model (number of elements and number of D.O.F.).

		80% 90% 100%						80% 90% 100%	0:17'32''	0:10'05''	0:6'49''	0:5'34''
	Time Average percentage	30% 40% 50% 60% 70%				Maxwell Transf. Interp. R-Adapt Preparation Remeshing	Time Sum percentage	30% 40% 50% 60% 70%	Maxwell Transf. Interp. R-Adapt Preparation Remeshing
		10% 20%	0:01'36''	0:01'17''	0:00'56''	0:00'48''		10% 20%
		0%						0%
			1	2	4	8			1	2	4	8
				Number of Procesors				Number of Procesors
				(a) Average

Table 4 .

 4 4 -Total times of the interpolation and transformation stages.

	Processors	1	2	4	8
	Forced	Int. Trans.	201.1(16%) 115.5(16%) 63.5(11%) 30.8(6%) 5.0(0%) 4.0(1%) 22.3(4%) 39.5(8%)
	Auto	Int. Trans	50.0(5%) 1.2(0%)	40.3(7%) 0.7(0%)	9.9(2%) 15.2(4%)	9.6(3%) 7.3(2%)

4.9. Study Cases purpose

  we have selected, as shown in fig.4.34, 3 point sensors and 3 elements at one of the symmetry faces of the model to read the output variables.

	Outer Element	Central Element
			Inner Element
	Outer Sensor
		Central Sensor	Inner Sensor

  • cut of the aluminum and copper tubes as shown in fig.4.43. For this simulation we assume symmetry at the sides cutting planes. Two cases are treated: The central tube is considered as

	(a) a rigid body
	(b) a deformable material

Table C .

 C 5 -Plastic coefficients for the Modified path-dependent Z.A. model. (Copper)[START_REF] Gould | A path-dependent constitutive model for gilding copper[END_REF]

	Parameter Value
	σ a [M P a]	50
	K[M P a]	0 (fcc materials)
	k[M P a]	1475
	S[K -1 ]	4.77e-5
	C 3 [K -1 ]	3.6e-4
	C 4 [K -1 ]	2.2e-5
	α[-]	0.88
	θ[M P a]	1048 + 10.42 ln( ε)
	ε0 [s -1 ]	2.1e11
	Σ[M P a]	289

Table C .

 C 6 -Plastic coefficients for the Takatsu model. (Aluminum 1050)[START_REF] Takatsu | High-speed forming of metal sheets by electromagnetic force[END_REF]

	Parameter Value
	K 1 [M P a]	118
	K 2 [M P a]	15.7
	n 1 [-]	0.27
	n 2 [-]	0.54
	ε0 [s -1 ]	1e-3

  RésuméLe formage électromagnétique est une technologie qui a gagné en intérêt dans les dernières décennies -grâce notamment à la formabilité accrue qu'il offre pour les matériaux à haute résistance et faible masse spécifique tels que les alliages d'aluminium et de magnésium. Un des défis majeurs au niveau du procédé réside dans la conception et l'étude au niveau de la pièce à fabriquer et sur l'interaction entre les différents aspects physiques: les ondes électromagnétiques comme source d'énergie, la thermomécanique contrôlant les évolutions de déformations et de contraintes, ainsi que l'étude de l'endommagement sous des sollicitations à grande vitesse. Ce travail est consacré à la mise au point d'un modèle et d'un outil numérique prédictifs capable de traiter l'interaction entre électromagnétisme et thermo-mécanique dans un cadre éléments finis en 3D. Nous introduisons les modèles de calcul pour l'électromagnétisme : l'approche permettant d'inclure la géométrie des pièces, le couplage électrique avec le générateur,

	entre autres. Nous poursuivons avec les
	techniques de calcul nécessaires pour
	coupler le calcul électromagnétique avec les
	calculs thermo-mécaniques -en mettant
	l'accent sur le problème du suivi des
	déplacements de la pièce déformable dans le
	module électromagnétique. Nous introduisons
	aussi certains aspects plus physiques du
	procédé	tels	que		les	phénomènes
	d'élimination du retour élastique ou encore
	l'adhésion des surfaces (soudage). Dans le
	dernier chapitre, nous présentons les
	installations expérimentales disponibles au
	laboratoire.	Une	méthodologie	pour
	l'identification des paramètres électriques
	définissant les machines et nécessaires pour
	effectuer la simulation est introduite.
	Mots Clés
	Formage magnétique,	formage	rapide,
	éléments finis, analyse d'images.

Compared to the slower forming technologies

Nuclear research division of General Dynamics Corporation

Explicit methods are largely used in fast dynamics because they will always provide a solution, but no guarantee is given on whether the solution found corresponds to the equilibrium state. Implicit methods on the other hand may have difficulties converging in the numerical

This is valid if we do not include thermal effects in the coil that could change the electric conductivity.

Meaning the velocities of the body relative to the tool

The filter is composed of: An attenuator factor at the oscilloscope CA = 10 19.94/20 and the filter of the probe CP = -1.968e8. The total coefficient is C = CA • CP = -3.9266e25.

It can change depending on the direction of propagation with respect to the latices of the crystal selected in the material

Acknowledgements

by AIRBUS HELICOPTER, BMAX, Avantis Engineering and our lab CEMEF. Is was because of their funding and support that I could carry out my research.

I have to thanks Francis FOURNIER, who helped me with the design and adaptation of our experimental machine.

Chapter 3. Computational Modelling of Electromagnetism

Test case

In the following we present a case that puts in evidence this distinction but more importantly, it tell us that for our edge-based problem of solving the electromagnetic equations, the in-house code outperforms the initial global assembly offered by the PETSc library.

The study case corresponds to a 2 • 3D cut of the EM forming machine available in our facilities. The geometry of the coil is described in fig. 3.20. It is a pancake-like coil with an internal diameter of 30mm and an external diameter of 132mm. A flat aluminium sheet is placed at 0.5mm distance from the top face of the coil. For the simulation, the two solids are surrounded by an air box of 60mm × 104mm. The boundary conditions are defined as follows: two magnetic symmetry planes at the cutting surfaces. The rest of the surfaces are considered free. A zero potential condition is prescribed at the output nodes of the inductor. For the input condition we consider an initial potential of V o = 2500V olts. The machine is characterized by a resistance R m = 40mΩ, inductance L m = 69nH and capacitance C m = 552µF . In order to apply the correct conditions to each turn of the coil we pre-multiply each face by the scaling factor

Compute relative coverage parameters fig. 4.21

Dpq -1 Compute Coverage ratio and resemblance factor

end end Rank old partitionings with respect to the new ones Initialization of the ranking variable: rank i = 1 for q ∈ [1, N partitions ] do if q > 1 then for j = q . . . 2, do if RF p q > RF p j-1 then rank j = rank j-1 rank j-1 = q end else rank j = q exit end end end end

Conclusions

In this chapter we have explored the global methodology needed in order to set up a coupled simulation between the electromagnetic module and the thermo-mechanical module. It was selected to work following a sequential coupling strategy that allows handling the physics in a individual fashion while allowing communication at each time-step.

The fully immersed approach used for the meshing of the electromagnetic problem introduced the need for mesh adaptation techniques enabling to track the displacement and deformation coming from the mechanical computations and impose it to the air mesh. The first technique consisted in an automatic meshing strategy. It combines a progressive-mesh size-by-boxes together with the skin depth phenomena in order to adapt the mesh size to the materials at hand. The main target of this approach is to concentrate most of the elements around the solids. Including the "Distance from coil" parameter into the meshing process allowed to reduce the area affected by the finer mesh size but as the part moves, and depending on the value assigned to this parameter, the affected zone is lost. We have evaluated the effect on the mesh size as we saw that it allows to reduce the number of unknowns of the system as the part gets farther from the coil. A question that has yet to be answered is the effect that it might have on the physical results and the accuracy. Never- The following step consists in performing the same optimization procedure but comparing the experimental results against the analytical model as well.

In this case we set as variables R eq and L eq . The work bench is composed of: the coil that has been previously defined, which is placed below an epoxy layer of 0.5mm thick in order to avoid contact with any other metallic object. A metal platform enables installation of different types of matrices for the forming process. There are currently three experiments that have been set-up as shown in fig. 5.7:

• Indirect forming, which uses an intermediate punch to transmit the kinetic energy.

• Free bulging experiment, where the plate is positioned directly over the coil and it is mechanically hold by a polymer matrix.

• Surface printing experiment to evaluate the degree of details that are imprinted on the material after impact A detailed description is given in what follows. Deformation under the given conditions can be evaluated by means of an image correlation software such as VIC3D. At the current moment we can only carry out a post-mortem analysis by comparison of images before and after the magnetic loading. However, this analysis gives access to a full map of the displacement deformation field and allows precise identification of the fractured zones.

In order to perform an accurate correlation of this experimental results with the simulation it would be needed to: In first place, use the full geometrical representation of the coil. The lack of perfect symmetry has a great influence on the sample. The fracture zone observed in fig. 5.11 is a direct consequence of this imperfection (the Lorentz force density in this area is most likely higher compared to other points on the same circumference). This leads to the second point, which is the need to include a damage model at high strain rates regimes.

Appendix A

Mathematical formulas and identities

A.1 Time derivative of the determinant of the Jacobian matrix, Euler formula

Demonstration: Given a transformation matrix

its determinant is given in inditial notation by

The time derivative of the determinant is computed as follows

It can be demonstrated that the term ∂ ∂x l J êl = 0.

Without loss of generality, lets isolate one of the partial derivatives. Take for

A.1. Time derivative of the determinant of the Jacobian matrix, Euler formula

End of the demonstration.

B.2 Lorentz Contraction from Hertz-Heaviside electrodynamics -Demonstration in the ω-k space

Suppose that the source of an EM field is moving with the velocity v in the vacuum. The field equations are

If the source is a point particle, the term ρ e can be eliminated by considering an arbitrarily large domain. By changing from the space-time domain to the wavenumber-frequency domain we get Appendix D

Shock wave mechanics

Most of the electromagnetic forming applications have an impact condition driving the process after the energy has been transmitted to the work piece. This means that a full comprehension of the process requires some basic knowledge on impact mechanics.

Shock wave mechanics was formally founded in the 19th century with the works of Rankine and Hugoniot. Rankine was the first to show that within the shock a non-adiabatic process must occur. Hugoniot showed that in the absence of viscosity and heat conduction, conservation of energy implies conservation of entropy in smooth regions and a jump in entropy across a shock [Salas 2007].

At the beginning the theory was mainly targeted to the study of shock propagation in gaseous or perfect elastic solids for the mathematical simplicity in defining the equations of state (EOS) relating pressure, temperature and density.

The physics of shock waves mechanics in condensed or solid matter saw its formal beginning in the 1950s with studies focused on finding the equations of state for solids in dynamic regimes [START_REF] Kanel | Shock-wave phenomena and the properties of condensed matter[END_REF]]. Research continues today with several studies in which nano-scale simulations together with laser pulsed shock experiments allows to follow in a detailed fashion the evolution and propagation of the shock wave within metal crystals, enabling the identification of the vast phenomena occurring within the material in such conditions.

In this appendix we mention a few concepts the we consider important in understanding the phenomena and that we can use for the study of the highspeed forming process.

Appendix E

Isoparametric Nedelec Element The tetrahedral represented in figure E.1 is characterized by

The interpolation functions are defined as in Jin 2002:

• Scalar linear interpolation functions

Abstract

Magnetic pulse forming is a technology that has gained interest in the last decadesthanks to the increased formability it offers for high-resistance-low weight ratio materials such as aluminum and magnesium alloys. One major complexity of the process lies in the design and study at the work piece level and the interaction between the several physical aspects involved: the electromagnetic waves as source of energy, the thermo-mechanics controlling the strain and stress evolution, as well as the study of fracture and damage under high-speed loading conditions. This work is dedicated to the development of a predictive model and computational tool able to deal with the interaction between the electromagnetism and the thermo-mechanics in a 3D finite elements frame work. We introduce the computational aspects of the electromagnetism, from the selected approach to include the geometry of the parts down to the coupling with the electric machinery behind the process. This is followed by the computational techniques needed to couple the electromagnetic computation to the thermo-mechanical one with a special focus on the problem of tracking the displacement of the deformable part within the electromagnetic module. We also introduce some aspects more related to the physics of the process such as the phenomena of elastic spring-back elimination and surface bonding (welding). In the last chapter we present the experimental facilities available at the laboratory. A methodology for identification of the electric parameters defining the machinery and needed to perform the simulation is introduced.