Diogo

Giles

Sylvie Lucas

Keywords: TIG welding, 3D cathode-arc-pool-anode model, multiphysics modelling, interface modelling, partitioned algorithms iii

his direction. Moreover, I thank Laurent for his strong support and help when things got difficult.

To my father Nadim and my mother Dany, how can I describe my love and gratefulness for your infinite support and confidence. Everything I

Résumé

Ce travail présente le développement et l'analyse des phénomènes physiques du procédé de soudage TIG (Tungsten Inert Gas), ainsi que les approches numériques nécessaires au couplage entre l'arc-plasma et le bain de soudage. L'objectif de l'étude est de développer un modèle 3D unifié pour la prédiction des variables clés de soudage en fonction des paramètres procédé de base.

Les développements mathématiques, implémentés numériquement via Cast εM (logiciel de simulation numérique développé par le CEA et basé sur la méthode des éléments finis) couplent un modèle de cathode, d'arc plasma et de bain de soudage de manière robuste. Effectivement, une dérivation formelle des conditions thermo-hydrauliques à l'interface fluide est effectuée de façon à intégrer à la fois le transfert de chaleur de domaines conjugués et le phénomène de changement de phase. Cette étape de dérivation permet de s'assurer que la cinématique, la dynamique et les termes thermiques sont proprement pris en compte au niveau des conditions à l'interface arc-bain. Ensuite, les modèles de l'arc-plasma et du bain de soudage sont numériquement couplés à travers leur interface. Pour coupler les modèles de conservation d'énergie, une approche algébrique de type variables mixtes est proposée et implémentée ce qui permet d'utiliser les variables adaptées à chaque domaine pour résoudre monolithiquement le système thermique. De plus, afin d'unir les modèles de conservation de masse et de quantité de mouvement des domaines de l'arc-plasma et du bain de soudage, un algorithme de Dirichlet-Neumann et un algorithme quasi-monolithique ont été implémentés. Leurs performances numériques sont analysées et discutées. Afin d'identifier les paramètres clés du modèle entièrement couplé, plusieurs études de sensibilité ont été effectuées. Au cours de ces études, l'influence de la géométrie de la cathode, du courant électrique en entrée, des hypothèses d'interfaces et de viscosité du bain sont discutées. Il ressort de ces études que l'arc-plasma est significativement lié aux singularités géométriques de la cathode, ce qui montre l'importance du choix de géométrie pour la mise en place de modèles couplés pour la simulation. Par ailleurs, l'influence du courant électrique et de la viscosité dynamique sur la thermo-hydraulique du bain mettent en évidence la prédominance des forces à l'interface arc-bain et de la force de Lorentz dans la formation du bain de soudage.

Finalement, une étude de vérification du modèle 3D est proposée, étudiée et analysée. Le modèle vérifié est ensuite utilisé pour simuler une configuration 3D entièrement couplée qui prend en compte les effets de déplacement. Une étape de comparaison expérimentale est également effectuée en utilisant des résultats expérimentaux issus de la littérature. Ce travail constitue donc une base solide pour la mise en place, dans un futur proche, d'un modèle 3D couplé exploitable pour le CEA.

M ots -clés : Soudage TIG, Modèle 3D cathode-arc-bain-anode, modélisation multiphysique, modélisation des interfaces, algorithmes partitionnés v To my parents Nadim and Dany, without whom I wouldn't have gotten this far.

To my siblings Marianne, Steven and Joyce, who have always supported me and believed in me.

To the family sweetheart, who brought us light when all was somber, to Tamara.

In loving memory of my dear friend Jassem.

Objective of this chapter

This chapter serves to introduce the reader to the problem statement of this thesis and to motivate and justify the approaches used throughout this work. To this end, the long term goal of the French Atomic Energy Commission 1 (CEA) concerning TIG welding is presented and is followed by a brief history of multiphysics modelling of TIG welding at CEA. Next, I present scaling arguments that motivate the physical and mathematical strategies employed in this work. Finally, to set the stage for the simulation of more general 3D TIG welding configurations, the reader is introduced to the sensitivity, verification and 3D weld displacement studies performed in this thesis.

Industrial context

TIG arc welding is an ubiquitous technique used to join metallic pieces together. An arc is used to generate heat at high temperatures, allowing for the melting of the welding zone of the workpiece. In order to understand the mechanisms at play in TIG welding, the constituents of the welding system are first identified. The system is grossly made up of: an electric cathode, an arc and an electric anode. These three components serve as the electric circuit that heat the arc and turn it into a high temperature heat source. To ensure a stable and controllable electrification of the arc, the cathode is placed in a nozzle that transports a shielding gas, serving as both the electrical medium to the arc and as a protective environment (against oxidation) for the molten metal in the weld pool. The cathode is generally made of doped tungsten (with thorium or lanthanium), the arc is made of an inert gas (argon or/and helium) and the anode is the target workpiece to be welded [START_REF] Leonard | Welding handbook: welding processes[END_REF]. Figure 1.1 schematises the above description of a typical TIG welding system. As for the technique itself, it is of particular interest to the nuclear industry due to the clean and 1 Commissariat à l'énergie atomique 1.2. INDUSTRIAL CONTEXT high quality welds it produces. With CEA being one of the principal nuclear research institutes in France, the commission strongly invests in the efforts to master the TIG welding technique. In fact, the long term goal of CEA is to develop a fully predictive and complete computer model that is capable of using typical welding parameters (arc height, inlet current, materials etc...) to calculate the thermo-mechanical responses and metallurgical characteristics of the welded workpiece. Working towards this goal, CEA currently invests in the modelling of the following: thermo-mechanical, thermo-metallurgical and magnetothermo-hydrodynamic approaches. While the latter is the subject of this thesis, the reader is referred to the work of Pommier [START_REF] Pommier | Stress relaxation cracking in AISI 316L-type austenitic stainless steels[END_REF] for his work on the thermo-mechanical phenomena of welding and to Baumard [5,6] for her work on the thermo-metallurgy of 316L steel. The magneto-thermo-hydrodynamic approach (also called the multiphysics approach) is based on the modelling of the electrodynamically forced TIG arc, the interaction of the arc with the liquid weld pool that forms in the workpiece and the resulting geometric form of the weld pool. Fundamentally, this means that the efforts put into multiphysics modelling of TIG welding are to master the predictibility of the generated thermal field of the workpiece when welding. Indeed, with a predictive tool capable of accurately reproducing the thermal field from basic welding parameters, the use of experimentally calibrated thermal sources in mechanical and metallurgical studies can eventually be avoided [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Pommier | Stress relaxation cracking in AISI 316L-type austenitic stainless steels[END_REF][START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF][START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF]5,6].

A brief history of multiphysics modelling of TIG welding at CEA

The primary focus of this thesis is to advance the multiphysics approach and is structured on the historic developments of CEA in this field. In fact, the general objective of this thesis is based on the state of the multiphysics models at CEA in the year 2018. Therefore, to contextualise the work presented in this thesis, the historical advancements on the multiphysics modelling of TIG welding is briefly discussed. The different works are listed in a chronological manner.

From 2005 to 2009:

The first multiphysics model at CEA was the result of Brochard's thesis [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF]. The work extensively covers the literature on the magneto-thermo-hydrodynamic phenomena used to model the TIG system. Brochard implemented a computationally efficient approach to model the electro-thermal phenomena between the cathode, the arc and the anode. Furthermore, his work proved, by using scaling arguments and asymptotic expansions, that the TIG arc should be modelled as a dilatational and low Mach number fluid. The developed TIG arc model was then coupled to the weld pool via a fixed, non-deformable interface. The developed 2D axisymmetric numerical model was verified with carefully chosen analytical solutions and the phyical model was validated by comparing simulation results to the literature and to in house experiments. The work of Brochard set up the first coupled 2D axisymmetric TIG model at CEA. The model was applicable to the study of spot configurations for weak arcs (where the arc-pool interface negligibly deforms), where the torch and workpiece remain stationary w.r.t the laboratory frame of reference. This served as the first step in the construction of a numerical tool capable of predicting the thermal field in the workpiece, by the use of typical welding parameters.

From 2009 to 2012:

To push the model further, Kong's thesis [START_REF] Kong | Modelisation 3D d'ecoulements avec surface libre pour le soudage a l'arc TIG[END_REF] investigated the deformation of the arc-pool interface. Kong's work focused on the implementation of a fluid model with a deformable free surface using a mesh tracking method. The purpose of the work was to implement and experimentally validate a free surface model driven by imposed pressure and surface tension phenomena. This set the stage for the incorporation of surface deformation effects in the modelling of TIG weld pools. The work of Kong added to the weld pool model the capability to capture surface deformation effects. This served to extend the range of physical phenomena added to the TIG welding model.

From 2012 to 2015:

Advancing the model further, Nguyen's thesis investigated the unification of all the dominant weld pool phenomena [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF]. Nguyen's work focused on the modelling of a 3D weld pool and is based on a decoupled approach where the TIG arc is replaced by thermal and mechanical boundary conditions. Thus, only the weld pool is incorporated into the computational domain. The work implements the enthalpy approach to resolve the phase change problem and focuses on the implementation of an efficient 3D electrodynamic solver for the weld pool. Furthermore, the work quantifies the impact that the deformable pool surface has onto the fluid dynamics of the system and thus onto the geometry of the weld pool. Nguyen's developed 3D weld pool model, that emulates the arc as a heat source using calibrated thermal and mechanical boundary conditions, is applicable to welding configurations for arcs that are weakly affected by weld pool dynamics.

Furthermore, the model is capable of efficiently capturing the thermal phase change effects in the pool.

Summarising the state of multiphysics modelling of TIG welding at CEA before the start of my thesis, it comprised of the following: 1) a 2D axisymmetric model that couples the cathode, arc and workpiece domains for fixed arc-pool interfaces; 2) a decoupled 3D weld pool model for deformable weld pool surfaces.

Problem statement and principal objectives

To approach the goal of CEA concerning TIG welding, a fully coupled model that incorporates all dominant welding phenomena is necessary. Furthermore, a predictive TIG welding model requires that 3D effects be captured by the fully coupled model because industrial welding invokes 3D thermal phenomena. Thus, to push the predictability of the TIG modelling efforts at CEA further, the coupling of the functionalities of the previous models and the extension of the coupled model to 3D is imperative. Therefore, the principal objective of this work is to identify, develop and implement coupling techniques to join the previous arcplasma and weld pool models, and ensure the mathematical consistency of 3D extensions to the coupled model. The developments are necessary to ensure that a wide variety of welding configurations can be studied. To facilitate the execution of the principal objective, it is arranged into the following goals:

1) an identification of both the dominant physical terms at the arc-pool interface and the thermal models best adapted to capturing the range of thermal phenomena in the TIG system.

2) Identification, development (when needed) and implementation of numerical methods best adapted to resolving the fully coupled thermo-hydraulic system is necessary.

3) To better comprehend the behaviour of the fully coupled model, the physical and numerical sensitivities of the model are to be identified and analysed. This is important if the dominant phenomena in 3D configurations are to be understood.

4) The verification of the the fully coupled 3D model is necessary for the model to be used in a reliable manner. 5) A 3D welding configuration with displacement effects is to be constructed, simulated and compared to an experimental case. Thus, to tackle the principal objective of this thesis I begin by presenting the model coupling methodologies that are used in this work.

Coupling the arc-plasma and pool: a scaling argument

Coupling the arc-plasma and the weld pool at their interface extends the applicability of the global model to more general welding configurations, where the arc-pool system can strongly interact. However, because the arc-plasma and the weld pool models are driven by different phenomena, the numerical implementation of the coupling requires attention. Furthermore, implementing the coupling of the arc and pool subdomains numerically should be compatible with the dominant thermo-hydraulics of the two subdomains 2 . In this thesis, compatibility is defined by the stability of the numerical algorithms and their physical coherence (i.e the algorithm should respect the conservation laws). To this end, identifying and analysing both the dominant dynamics at the interface and the adapted thermal models between the arc and pool subdomains is paramount. Thus, in looking to couple the two models, scaling arguments are requisite, and the scaling methodology follows that described by Ruzicka [START_REF] Mc Ruzicka | On dimensionless numbers[END_REF]. Before presenting the scaling arguments, the governing arc-pool momentum and energy conservation equations at the fluid interface are invoked. The reader is referred to sections 2.4 and 2.5 for a detailed presentation and discussion of the conservation laws. where n, W e, P , Bo, ∆ρ, Ca et σ represent the interface normal vector, Weber number, fluid pressure, Bond number, difference in density, Capillary number and the normal viscous stresses, respectively. The non-dimensional numbers, as they apply at a typical arc-pool interface, are described and expressed in table 1.1. Furthermore, figure 1.2 schematises the interaction of the arc and pool at their interface for 2 The coupling of the mesoscopic electrodynamic models of the two subdomains is straightforward and thus is not considered in the analysis. an arbitrary TIG spot configuration, which helps simplify the following scaling arguments 3 . The main stagnation zones represented in the figure as a, b and c refer to zones of stagnating arc or pool flow where the stagnation pressure of the fluid rises as a consequence. Indeed, because the arc-pool system is a two fluid structure with dynamics along different spatial scales, the three presented stagnation zones manifest stagnation pressures of different magnitude. In fact, stagnation zone a is dictated by the arc flow, while zone b is dictated by the fast moving interface velocites and zone c by the slower resurfacing flow in the pool. Moreover, according to elementary fluid dynamics theory, the stagnation pressure can be scaled to P ∼ ρU 2 [1], which translates to the dependence of the Weber number onto the square of the characteristic velocity of the zone of interest W e = ρU 2 /(γL). Centering the scaling arguments in the vicinity of the stagnation zones, the following characteristic scales and quantities are chosen as they manifest for a typical welding configuration [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF][START_REF] Mills | Factors affecting variable weld penetration[END_REF][START_REF] Morohoshi | Surface tension of liquid iron as functions of oxygen activity and temperature[END_REF]:

Arc              U a = 100 m • s -1 µ a = 10 -4 Pa • s L a = 10 -3 m ρ a = 10 -1 kg • m -3 Pool                                        U p = 5 × 10 -1 zone b 10 -2 zone c m • s -1 µ p = 10 -3 Pa • s L p = 10 -2 m H z = 10 -3 m ∆T = 500 K ρ p = 7000 kg • m -3 | dγ dT | = 10 -4 N • m -1 • K -1 γ = 2 N • m -1
where the pool characteristic velocity U p has two dominant scales that interact with the arc-pool interface at the stagnation zones b and c. These stagnation zones significantly impact the pool Weber number, due taken into account when coupling the arc and pool momentum models because the characteristic SAP and M P S numbers are of similarly significant magnitudes. As for the dynamics along the normal of the arc-pool interface, the arc and pool pressure need to be similarly considered. The Weber number of the arc being relatively significant at stagnation zone a w.r.t the pool Weber number at zone c implies that the arc dominates the deformation of the interface downstream the arc axis. Furthermore, the relatively significant pool Weber number at zone b indicates that the pool pressure dominates any potential humping of the interface. The Bond and arc Capillary number being relatively equivalent in their contribution to interface deformation, implies the importance of incorporating both the gravity and arc viscous stress forces. Lastly, regardless of the miniscule magnitude of the pool Capillary number, the viscous contribution of the pool is included into the coupling effort of the arc-pool models so as to respect viscous flow theory [1].

Next, a simple dimensionless thermal model is presented where only the ratio of internal energy contributions is discussed, and is expressed as:

( dT dt ) * + 1 Ste df dt = 1 P e (∇ 2 T ) * (1.3)
where Ste = c p ∆T /L f is the Stefan number, f the phase fraction and P e the Peclet number. The Peclet number is significant in both the arc and pool models as was shown in Brochard's work [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF], which justifies the incorporation of heat transport effects in both subdomains. Moreover, the Stefan number is the ratio of sensible heat c p ∆T to the latent heat of fusion L f , and is only valid if the phase fraction, f cnst is not a constant. The phase fraction represents the percentage of one phase over the other in a phase change problem (see equation (2.26) in chapter 2). Thus, in looking to couple the thermal model of the arc-plasma to that of the weld pool, the characteristic Stefan number of both subdomains is calculated using characteristic weld pool specific and latent heat values (see appendix A).

Arc f = 1 implies Ste is undefined Pool 0 ⩽ f ⩽ 1 -→ Ste = c p ∆T L f ≈ 600 × 500 2.5 × 10 5 ≈ 1
The Stefan number being undefined in the arc means that the implemented arc thermal model is not required to be adapted for phase change effects. Choosing not to model metal evaporation effects in the arc-plasma means that the temperature based model of Brochard, which is adapted to capturing the electro-thermal effects [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF], should be used. However, with Ste ≈ 1 in the weld pool, the sensible and latent heat contributions are of the same order, which requires that the implemented thermal model be adapted to capturing phase change effects. To ensure that both the sensible and latent heat effects are captured implicitly by the pool thermal model, an enthalpy based thermal model is used [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF][START_REF] Vaughan R Voller | An enthalpy method for convection/diffusion phase change[END_REF][START_REF] Nedjar | An enthalpy-based finite element method for nonlinear heat problems involving phase change[END_REF].

Having identified the dominant dynamic forces at the arc-pool interface and the type of thermal models necessary to capture both electro-thermal and phase change effects, the identification of the adapted numerical coupling techniques is necessary. The coupling algorithms implemented in this thesis are argued for in the following:

Mass and momentum coupling algorithm:

Recalling that the absolute magnitudes of the SAP and M P S numbers are significant as opposed to the weaker absolute magnitudes of the Weber, Bond and Capillary numbers, I base the structure of the coupling algorithm on the physical characteristics of the interface. Thus, the algorithm should ensure that a convergent and conservative communication of the interfacial forces is respected while allowing for a mass conserving deformation at the interface. To this end, numerically coupling the mass and momentum equations between the two subdomains is tested using two methods, a partitioned Dirichlet-Neumann and a quasi-monolithic algorithm. The partitioned Dirichlet-Neumann algorithm transmits, at every iteration of a fixed mesh, the arc pressure and shear forces to the pool as Neumann boundary conditions, and the pool surface velocities to the arc as Dirichlet boundary conditions. This procedure is iterated until force/velocity convergence at the fixed interface is achieved. As for the quasi-monolithic algorithm, it ensures a conservative communication of the fixed interface forces without requiring internal iterations.

The algorithm uses the Lagrange Multiplier method [START_REF] Pellet | Dualisation des conditions aux limites[END_REF] to ensure the conservation of forces at the fixed interface. After either algorithm is executed on a fixed mesh, the normal interfacial forces and a mass conservation constraint are used to calculate the deformation at the interface. Once the interface deformation field is calculated, a mesh deformation algorithm based on Huang's work [START_REF] Huang | Mathematical principles of anisotropic mesh adaptation[END_REF] is used. The mesh deformation method conserves mass and allows for the capture of the weak deformation of the arc-pool interface. This approach avoids the use of interface tracking methods which in this thesis are considered excessive. For example, the volume of fluid, phase field and level set methods are computationally costly and can introduce issues with surface tension and mass conservation [START_REF] Cadiou | Modélisation magnéto-thermohydraulique de procédés de fabrication additive arc-fil (WAAM)[END_REF]. The use of interface tracking methods are thus reserved for problems where violent interface deformation phenomena (for W e >> 1

and Ca >> 1) such as interface detachement/reattachement and when drop formation phenomena occur [START_REF] Cadiou | Modélisation magnéto-thermohydraulique de procédés de fabrication additive arc-fil (WAAM)[END_REF][START_REF] James | Level set methods for fluid interfaces[END_REF]. The mass and momentum coupling algorithms are thus developed in section 3.6.

Conjugate heat transfer coupling algorithm:

Recalling that the Stefan number is undefined in the arc-plasma but is Ste ≈ 1 in the weld pool, I base the structure of the coupling algorithm on the thermal characteristics of the two subdomains. In order to both respect energy conservation and be compatible with conjugate heat transfer and phase change effects in the arc-pool system, a mixed variable monolithic algorithm is proposed and implemented. The mixed variable algorithm couples the temperature based heat transfer model of the arc-plasma [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF] to the enthalpy based heat transfer model of the weld pool [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF]. This allows the algorithm to simultaneously capture both the electro-thermal and heat transport phenomena of the arc, and the phase change phenomena present in the weld pool. The algorithm uses the Lagrange multiplier method [START_REF] Pellet | Dualisation des conditions aux limites[END_REF] to communicate the heat flux across the arc-pool interface, by weakly imposing the temperature continuity using the temperature field of the arc and the enthalpy field of the pool. The monolithic mixed variable algorithm strongly conserves energy while also satisfying temperature continuity across the interface within a couple of global iterations. These methods are developed in section 3.5.

Having presented the choice of the coupling algorithms that are compatible with the dominant dynamics and thermal phenomena between the arc and pool subdomains, I briefly note the range of applicability of the coupled model. The studied unified model in this work is stationary in that the global algorithm searches for steady state solutions directly without integrating in time. This choice is taken because it is imperative that the coupling algorithms are first developed to respect the dominant conservation laws of the physical problem and converge numerically [START_REF] Viot | Couplage et synchronisation de modeles dans un code scenario d'accidents graves dans les reacteurs nucleaires[END_REF][START_REF] Patricot | Couplages multi-physiques : évaluation des impacts méthodologiques lors de simulations de couplages neutronique/thermique/mécanique[END_REF][START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF][START_REF] Carlos A Felippa | Partitioned analysis of coupled mechanical systems[END_REF]. Thus, I focus on ensuring a robust coupling of the arc-pool system and reserve the implementation of a time integrator for a future work. Furthermore, the incorporation of filler metal into the arc-pool system is also reserved for a future work.

Thus, to continue tackling the principal objective of this thesis I briefly motivate the sensitivity studies and the model extension to 3D that is performed in this work.

An extension to 3D configurations

Once a robust and conservative coupling algorithm that includes the dominant phenomena in the TIG welding system is implemented, the model can be generalised and extended to 3D configurations. However, before extending the model to 3D, identifying and understanding the sensitivity of the unified model to the choice of the coupling algorithm and to certain geometric and physical parameters is necessary.

This is important if the physical and numerical implications of the proposed model are to be understood before setting up complex 3D configurations. To this end, five sensitivity studies are performed on 2D axisymmetric configurations. Setting up the studies on axisymmetric configurations allows for good tradeoff between acceptable physical representativity and affordable computational costs. The sensitivity studies performed are listed and described in the following:

1) The sensitivity of the arc-plasma to the cathode geometry is studied. The cathode geometry is identified as an often overlooked essential parameter in the numerical simulation of TIG welding.

However, it is essential because the electron emission zone at a cathode can be the most influential boundary condition as it is the source of momentum for the arc [START_REF] Hsu | Study of the free-burning high-intensity argon arc[END_REF][START_REF] Delgado-Álvarez | Dimensionless representation of the column characteristics and weld pool interactions for a DC argon arc[END_REF]. A TIG arc is thus sensitive to the variations in cathode shape [START_REF] Leonard | Welding handbook: welding processes[END_REF]. The sensitivity of these arcs to cathode truncation angle has been studied both experimentally and numerically by multiple authors [START_REF] Goodarzi | The effect of the cathode tip angle on the GTAW arc and weld pool: I. mathematical model of the arc[END_REF]32,[START_REF] Alber | Gas-tungsten-arc cathode and related phenomena[END_REF][START_REF] Matsuda | Study on gas-tungsten-arc electrode (report 1): Comparative study of characteristics of oxide-tungsten cathode (welding physics, process & instrument)[END_REF][START_REF] Niles | Weld thermal efficiency of the GTAW process[END_REF][START_REF] Mills | Factors affecting variable weld penetration[END_REF]. Moreover, disagreement between numerical and experimental results (see Lago, Goodarzi, Tsai, Sadek [START_REF] Lago | A numerical modelling of an electric arc and its interaction with the anode: Part I. the two-dimensional model[END_REF][START_REF] Goodarzi | The effect of the cathode tip angle on the GTAW arc and weld pool: I. mathematical model of the arc[END_REF][START_REF] Tsai | Heat distribution and weld bead geometry in arc welding[END_REF][START_REF] Alber | Gas-tungsten-arc cathode and related phenomena[END_REF]) motivate a more careful analysis of the effect of cathode geometry onto TIG arcs. For example, Tsai [START_REF] Tsai | Heat distribution and weld bead geometry in arc welding[END_REF] shows, experimentally, that convective heat transfer to the work piece increases as cathode truncation angle decreases, and that is not evident in the numerical results of Goodarzi [START_REF] Goodarzi | The effect of the cathode tip angle on the GTAW arc and weld pool: I. mathematical model of the arc[END_REF]. Furthermore, not many numerical investigations into the effect of cathode tip size and shape, which can have strong effects on the arc [START_REF] Alber | Gas-tungsten-arc cathode and related phenomena[END_REF], have been performed, as authors generally choose arbitrary pointed and chamfered tip shapes (see [START_REF] Goodarzi | The effect of the cathode tip angle on the GTAW arc and weld pool: I. mathematical model of the arc[END_REF][START_REF] Bini | Numerical and experimental study of transferred arcs in argon[END_REF]); thus motivating this work. This study is presented in chapter 4.

2) The sensitivity of the model to the two implemented mass and momentum coupling algorithms is studied. The Dirichlet-Neumann and the quasi-monolithic algorithms are tested for three different TIG spot configurations, where the convergence rate, solution sensitivity and total iteration count are compared and analysed. This algorithm sensitivity study helps identify the algorithm that is better adapted to resolving the unified TIG model. This serves to ensure that the chosen algorithm is robust and relatively computationally efficient, which will facilitate computational efforts for 3D

configurations. This study is presented in section 5.6.1.

3) The sensitivity of the weld pool to the inlet electric current is studied. Although the influence of inlet electric current onto weld pools is experimentally understood and quantified (see Mills et al. [START_REF] Mills | Factors affecting variable weld penetration[END_REF] for an extensive review), the sensitivity of the driving forces that dictate the pool dynamics is not well established [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Mills | Factors affecting variable weld penetration[END_REF]. In general it is known that the increase of inlet electric current increases both the available heat and flow speeds in the arc which increases pool width and depth [START_REF] Mills | Factors affecting variable weld penetration[END_REF][START_REF] Tsai | Heat distribution and weld bead geometry in arc welding[END_REF][START_REF] Leonard | Welding handbook: welding processes[END_REF].

However, the sensitivity of the interfacial forces has not been thoroughly discussed in the literature and this has lead to certain authors omitting either the interfacial shear or normal force contributions [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF][START_REF] Kong | Modelisation 3D d'ecoulements avec surface libre pour le soudage a l'arc TIG[END_REF][START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF]. In fact, the similar magnitudes of the SAP and M P S (see section 1.3) numbers 1.4. AN EXTENSION TO 3D CONFIGURATIONS at a typical arc-pool interface indicate the competing roles the shear and Marangoni forces have in driving the dynamics of the weld pool. Furthermore, the use of augmented liquid steel viscosity values by some authors [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF][START_REF] Traidia | Multiphysics modelling and numerical simulation of GTA weld pools[END_REF][START_REF] Kim | Heat and fluid flow in pulsed current gta weld pool[END_REF] dampens the influence of the inlet electric current onto the pool geometry. Thus, a discussion of the sensitivity of the leading driving forces, dependent on inlet current, in the pool offers insight that allows for clearer interpretations of the dynamics of more complex welding configurations. This study is presented in subsection 5.6.2.1.

4) The sensitivity of the pool system to the chosen interface hypothesis: deformable or fixed/nondeformable, is studied. In fact, the low magnitudes of the Weber numbers at the arc-pool interface (see section 1.3) have been used to justify the use of fixed interfaces in the welding models [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Kidess | The influence of surfactants on thermocapillary flow instabilities in low prandtl melting pools[END_REF].

However, the work on weld pool dynamics in Nguyen's and Traidia's PhDs shows the significant influence the pool surface has onto the pool geometry [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF][START_REF] Traidia | Multiphysics modelling and numerical simulation of GTA weld pools[END_REF]. Thus, a discussion of the sensitivity of the pool geometry to the interface hypothesis offers insight onto and allows for clearer interpretations of the dynamics that influence pool penetration. This study is presented in subsection 5.6.2.

5)

The sensitivity of the arc-pool system to different pool viscosities is studied. Experimentally, the viscosity of liquid 316L steel is consistently reported by multiple authors [START_REF] Bobkov | Thermophysical properties of materials for nuclear engineering: a tutorial and collection of data[END_REF][START_REF] Choong | Thermophysical properties of stainless steels[END_REF][START_REF] Morville | Modélisation multiphysique du procédé de Fabrication Rapide par Projection Laser en vue d'améliorer l'état de surface final[END_REF][START_REF] Kidess | The influence of surfactants on thermocapillary flow instabilities in low prandtl melting pools[END_REF]; however, much of the literature on numerical simulations of weld pool dynamics uses inconsistent values (see [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF][START_REF] Kong | Modelisation 3D d'ecoulements avec surface libre pour le soudage a l'arc TIG[END_REF][START_REF] Traidia | Multiphysics modelling and numerical simulation of GTA weld pools[END_REF][START_REF] Yau | Modélisation numérique instationnaire pour la simulation du soudage TIG avec couplage plasma/bain de fusion[END_REF][START_REF] Kim | Heat and fluid flow in pulsed current gta weld pool[END_REF]). In fact, the differences in the experimentally reported and numerically implemented pool viscosity values are 10-fold in some instances. Thus, it is paramount that the influence of pool viscosity onto pool geometry and arc-pool dynamics be studied and analysed. This study is presented in subsection 5.6.2.3.

Once the sensitivities in the unified TIG model are analysed and understood, the extension of the model to 3D configurations can be performed. To extend the unified model in a coherent manner the mathematical consistency and physical validity of the model must be argued for. To this end, the verification and validation of the extended 3D model is performed. The developments of the 3D model are listed and described in the following:

1) The verification study consists of comparing the simulation results of a 2D axisymmetric TIG spot configuration to that of its 3D Cartesian analogue. This serves to verify the mathematical model and serves to indicate whether any mathematical inconsistencies have been programmed into the extended model. First, the rotational symmetry of the fields calculated for the 3D spot configuration is verified. The verification of the rotational invariance of the fields is necessary so as to be able to use the 2D axisymmetric case as a reference, and to compare the comparable. Next, the rotationally invariant 3D fields are compared to the 2D axisymmetric reference simulation and the the spatial convergence of the 3D fields is briefly discussed. This study is presented in section 6.1.

2) Lastly, the verified 3D model is compared to an experimental case for preliminary validation. The study includes weld displacement effects and is compared to an experimental configuration taken from Koudadje's work [START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF]. Effectively, the displacement effects are modelled by assigning to the workpiece a displacement velocity w.r.t the laboratory observer, while the cathode is held fixed. The simulation results of the fully coupled 3D cathode-arc-workpiece weld displacement model are then discussed and analysed. Furthermore, the calculated weld pool is compared to experimental macrographic images in terms of geometry and size. This study is presented in section 6.2.

Thus, with insight about the sensitivities of the unified model to multiple physical and numerical parameters and with a verified and generally valid 3D model, the model can be used to study different 3D welding configurations.

Introduction

In this chapter the mathematical models used to enforce the dominant electrodynamic laws, mass, momentum and energy conservation laws in the TIG cathode-arc-workpiece system are presented. The chosen arc and workpiece models are compatible with both conjugate heat transfer and phase change phenomena, while also allowing for the capture of the deformable interface phenomenon between the arc and the weld pool. Furthermore, the thermo-hydraulic equations at the arc-pool interface are formally derived, bringing into context the importance of non-negligible advective terms that could govern the interface thermo-hydraulics in some cases.

The chapter begins with a description of the relevant geometric domains in section 2.2, following with a presentation of the mesoscopic electrodynamic model in 2.3, the modelled momentum and mass conservation laws in section 2.4 and the modelled energy conservation laws in section 2.5.

Geometric preliminaries

The geometric construction of the TIG system, as discussed below, is schematised in the figure 2.1. To facilitate the mathematical discussion performed in this chapter, key geometric elements are presented.

Representing the global TIG welding system Ω T IG mathematically:

Ω T IG = Ω cat ∪ Ω pla ∪ Ω ano
where Ω cat is the cathode subdomain, Ω pla the arc plasma subdomain and Ω ano the anode/workpiece subdomain. The subdomain intersections, being important for the assignment of interfacial conditions, are identified as:

Γ CP I = Ω cat ∩ Ω pla Γ AP I = Ω pla ∩ Ω ano
where the subdomains of intersection, Γ CP I and Γ AP I , are aptly named the cathode-plasma interface (CPI) and the anode-plasma interface (API), respectively. The Γ CP I is a rigid, undeformable subdomain that consists of an intersection between a solid and a fluid; while the Γ AP I , consists of both a rigid, undeformable zone as well as a deformable zone1 . The Γ AP I connects the intersection of the arc-plasma to both the solid and the liquid weld pool sections of the workpiece. Furthermore, the boundary of the global TIG welding system Ω T IG can be defined, and is mathematically represented as:

∂Ω T IG = ∂Ω cat ∪ ∂Ω pla ∪ ∂Ω ano
but for which:

∂Ω ano ∩ ∂Ω pool = ∅
where the weld pools studied in this thesis are not deep enough to perforate the anode nor large enough to intersect with its lateral boundaries.

Electrodynamic equations

The TIG welding system is made up of a closed electric circuit made up of a cathode, a resistor and an anode. The cathode is the source of both the electric current and electric potential, while the resistor transforms input electric power into high quality heat2 and simultaneously transports the heat to the anode. The resistor in the TIG circuit is the arc-plasma, while the anode is the workpiece. In fact, it is because the electric circuit both generates and transports heat to the workpiece that the electrodynamic phenomena of the TIG system are the main drivers of the welding process. While the arc-plasma and more specifically its interaction with the electrode boundaries are fundamentally governed by complex multi-specie electrodynamics [START_REF] Lowke | LTE diffusion approximation for arc calculations[END_REF]2,[START_REF] Gleizes | Thermal plasma modelling[END_REF] the main driving phenomena can be captured by modelling the major mesoscopic phenomena [START_REF] Tanaka | A unified numerical modeling of stationary tungsten-inert-gas welding process[END_REF][START_REF] Tanaka | An introduction to physical phenomena in arc welding processes[END_REF][START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Gleizes | Thermal plasma modelling[END_REF][START_REF] Lago | A numerical modelling of an electric arc and its interaction with the anode: Part I. the two-dimensional model[END_REF]. The mesoscopic electrodynamic phenomena respects both the local thermodynamic equilibrium (LTE) hypothesis and the continuum model for the modelled solids, fluids and their interfaces [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF]. However, because the arc-plasma at the electrode boundaries significantly deviates from LTE as a result of its cooler boundary temperatures, modified electrical conductivities are required for those regions. To account for the non-LTE electrical conductivity in the electrode boundary regions the LTE diffusion approximation method from Lowke et al. is used [START_REF] Lowke | LTE diffusion approximation for arc calculations[END_REF]. This method assigns, in the cathode region, a modified conductivity taken at the value of the arc-plasma at a δl cat = 0.1 mm Euclidean distance from the CPI; while in the anode region, the modified conductivity value is assigned to the value at a δl ano = 0.4 mm Euclidean distance from the API (see section 3.3.1 for details). The mesoscopic electrodynamic phenomena are thus modelled by an electrostatic and magnetostatic approach adapted from Brochard's work [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF].

Electrostatic model

The equation used to resolve the electrostatic model for the electric potential ϕ across the entire domain is as follows:

-∇ • σ * ∇ϕ = 0 (2.1)
where the electrical conductivity σ * is represented by the following piece-wise function:

σ * (T ) =        σ * cat (T ) r ∈ Ω cat σ * pla (T ) r ∈ Ω pla σ * ano (T ) r ∈ Ω ano
where the electrical conductivities functions of both temperature and material type are given in appendix A.

Magnetostatic model

The equations used to solve the magnetostatic model for the magnetic field B across the entire domain are as follows:

∇ ∧ B = µ 0 j (2.2) ∇• B = 0 (2.3) for j = -σ * ∇ϕ (2.4)
where µ 0 = 4π × 10 -7 N•A -2 is the magnetic permeability in vacuum, assumed constant in the entire system. The set of magnetic equations (2.2) and (2.3) are defined as Ampere's law and Gauss's law, respectively [START_REF] David | Classical electrodynamics[END_REF].

Electromagnetic sources

Both the heat generation and its transport to the anode in the TIG system are modelled by electromagnetic phenomena. The transformation of electrical power to heat in the TIG system is largely dominated by the Joule effect [START_REF] Tanaka | An introduction to physical phenomena in arc welding processes[END_REF] and is modelled as:

s Joule = σ * ∇ϕ • ∇ϕ (2.5)
where the phenomenon is present in all three subdomains, but is strongest in the arc-plasma. According to multiple authors [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Tanaka | A unified numerical modeling of stationary tungsten-inert-gas welding process[END_REF][START_REF] Lago | A numerical modelling of an electric arc and its interaction with the anode: Part I. the two-dimensional model[END_REF][START_REF] Gleizes | Thermal plasma modelling[END_REF], the transport of this heat from the hot arc-plasma to the colder anode is driven by the dominant source of arc momentum, the Lorentz force. Furthermore, this force also contributes as a source of momentum in the weld pool, although to a much lesser extent [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF]. The Lorentz force, f Lor , is modelled as:

f Lor = j ∧ B = -σ * ∇ϕ ∧ B (2.6)
and couples to the momentum equations in the arc-plasma and weld pool; as will be presented in the following sections.

Boundary and interface conditions

With the electrostatic and magnetostatic equations defined as three partial differential equations, the definition of appropriate closure relations defined at the dim(Ω) -1 subdomains is necessary. Indeed, if Ω ∈ R 2 then the closure relations would apply on all (∂Ω, Γ) ∈ R 1 ; and similarly, if Ω ∈ R 3 then the closure relations would apply on all (∂Ω, Γ) ∈ R 2 . Furthermore, because the electrodynamic model is chosen to obey the mesoscopic LTE hypothesis and because the magnetic permeability coefficient is constant in the entire domain, continuity for both the electric potential and magnetic field is maintained at the CPI and API. The generalised boundary and interface conditions for equations (2.1), (2.2) and (2.3) are presented.

Boundary conditions

The electrostatic boundary conditions are governed by an imposed inlet current at the cathode and an electric ground at the anode. The lateral boundaries, being both cold and far from the electric current source are considered as electrically insulated where current cannot pass through. The insulation condition also applies for axes or planes of symmetry, if they are defined. This is mathematically expressed as:

Inlet: -σ ∂ϕ ∂n ∂Ω in = j imp = I imp ∥∂Ω in ∥ • n (2.7) Ground: ϕ ∂Ω ground = 0 (2.8)
Insulated or symmetry:

-σ ∂ϕ ∂n ∂Ω insu • n = 0 (2.9)
where j imp and I imp are the imposed current density and the imposed current, respectively. As for the magnetic boundary conditions, Ampere's law from equation (2.2) requires only a parametrisation of the spatial dimension because it is fundamentally a circulation, and is dependent only on the electric current inside the domain of interest. As for Gauss's law from equation (2.3), it requires at least one boundary condition to fix the magnetic field [START_REF] David | Classical electrodynamics[END_REF][START_REF] Jiang | The least-squares finite element method: theory and applications in computational fluid dynamics and electromagnetics[END_REF]. Furthermore, for symmetry axes and planes, the tangential field to the plane is null. The TIG circuit is assumed isolated and sufficiently far from any external magnetic field, thus:

Fixed:

B ∂Ω f ixed = 0 (2.10) Symmetry: B • t ∂Ωsym = 0 (2.11)
where t is the total tangential vector and manifests in 3D as t = t 1 + t 2 . Furthermore, the magnetic symmetry condition in equation (2.11) is chosen to be perpendicular to the electric current symmetry condition in equation (2.9).

Interface conditions

The simplified electrodynamic model makes for continuous electric potential and magnetic fields across both the CPI and the API. This is mathematically expressed for the electric potential as: at the CPI:

ϕ Γcat = ϕ Γ pla (2.12)
at the API:

ϕ Γ pla = ϕ Γano (2.13)
while for the magnetic field this is expressed as: at the CPI:

B Γcat = B Γ pla (2.14)
at the API:

B Γ pla = B Γano (2.15)

Momentum and mass conservation laws

The TIG welding system is generally made up of two fluid subdomains, identified as the arc plasma and the liquid weld pool. The arc is the subdomain that transports energy from the cathode to the workpiece, while the weld pool redistributes its absorbed heat in the target weld zone. If both the heat transport phenomenon of the arc and the redistribution effect3 of the pool are to be captured by the implemented model, the flow dynamics of both subdomains are to be necessarily treated. To be able to predict the flow dynamics of the system, both the conservation of mass and momentum are necessary. Therefore, the equations chosen to model the fluid conservation laws are based on the dilatational fluid model of Brochard [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF] for the arc-plasma and the incompressible fluid model of Nguyen [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF] for the weld pool.

Like the dilatational model from Brochard [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF], the arc-plasma is not considered to dilate like an ideal gas, rather it dilates according to the experimentally determined density-temperature relation reported by

Boulos [START_REF] Maher I Boulos | Thermal plasmas: fundamentals and applications[END_REF]. The main hypotheses taken in modelling the two subdomains are discussed in the following subsections.

Dilatational arc-plasma fluid model

The arc-plasma is modelled as a fluid under the following hypotheses, which are argued for in detail in Brochard [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF] and considered valid for the entirety of this thesis:

Newtonian fluid: The arc-plasma has isotropic thermophysical properties, independent of the strain rate.

Dilatational fluid: The arc-plasma is considered to behave like a dilatational fluid in that the dependence of the fluid density onto its temperature cannot be neglected. This assumption is argued for because the argon gas that transforms into the arc-plasma downstream the cathode has temperature values that can range from room temperature4 to the tens of thousands of Kelvin in the arc core.

This hypothesis requires that both the density of the arc-plasma and its dynamic viscosity be strongly dependent onto temperature. Furthermore, it is assumed that any variations in the dynamic pressure of the fluid are asymptotically small when compared to the thermodynamic pressure of the fluid.

This decouples the dependence of the fluid density onto the fluid pressure. This hypothesis is argued for because the arc-plasma in TIG welding applications flows at low Mach numbers.

Laminar flow: For typical TIG welding configurations, the flow of the arc-plasma is assumed to remain laminar. This hypothesis is considered because the flow of the arc-plasma is at Reynolds numbers that do not exceed Re ≲ 10 3 .

Resistive magnetic fluid:

The arc-plasma is assumed to be magnetically resistive in that the magnetic field is decoupled from the fluid flow and no magnetic advection occurs. This hypothesis is considered due to the low magnetic Reynolds number of the arc-plasma Re m << 1.

Having discussed the main hypotheses taken in the fluid model of the arc-plasma, the modelled mass and momentum conservation laws are presented in the following.

Conservation of mass:

Before presenting the equation modelling the conservation of mass, the following presentation relates the dependence of the density of the arc-plasma to its temperature. The density of the arc-plasma is defined as a thermodynamic state variable, and is dependent on the pressure and temperature [START_REF] Theodore L Bergman | Fundamentals of heat and mass transfer[END_REF], in:

dρ = ∂ρ ∂P dP + ∂ρ ∂T dT (2.16)
where P is the arc-plasma pressure. However, for a dilatational fluid that assumes a negligible dependence of its density onto its pressure, equation (2.16) becomes:

if ρ = f (T ) → ∂ρ ∂P ≈ 0 ; dρ = ∂ρ ∂T dT (2.17)
where both the density and temperature state variables are defined over the arc-plasma subdomain (ρ, T ) ∈ Ω pla . This allows for the use of the following expansions:

dT = ∂T ∂x dx + ∂T ∂y dy + ∂T ∂z dz = ∇T • dr dρ = ∂ρ ∂x dx + ∂ρ ∂y dy + ∂ρ ∂z dz = ∇ρ • dr
where by combining the above expansions to equation (2.17), the following expression becomes valid:

∂ρ ∂T = dρ dT = ∇ρ • dr ∇T • dr → dρ dT = ∇ρ • 1 ∇T (2.18)
where the derivative of density w.r.t temperature is defined across the entire Ω pla subdomain. To relate equation (2.18) to the flow of the arc-plasma, the law of conservation of mass for fluids is introduced, in:

∇ • (ρu) = ρ∇ • u + ∇ρ • u = 0 → ∇ • u = - 1 ρ ∇ρ • u (2.19)
where u ∈ Ω pla is defined as the velocity field in the arc-plasma. Combining equations (2.18) and (2.19), the following equation is obtained:

∇ • u = - 1 ρ dρ dT ∇T • u = - d ln(ρ) dT ∇T • u (2.20)
where the d ln(ρ) dT is calculated from the density-temperature relation of argon at atmospheric pressure, as reported by Boulos [START_REF] Maher I Boulos | Thermal plasmas: fundamentals and applications[END_REF], and as plotted in appendix A. Interpreting equation (2.20), any thermal dilatation of the fluid is balanced by an increase in the flow rate away from the zone of fluid expansion.

Conservation of momentum:

The equation used to model the conservation of momentum in the dilatational arc-plasma is as follows:

ρu • ∇u = -∇P + ∇•µ ∇u + ∇ t u - 2 3 ∇• uI + f Lor (2.

Incompressible weld pool fluid model

The weld pool is modelled as a fluid under the following hypotheses, which are used in Nguyen [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF] and considered valid for the entirety of this thesis:

Newtonian fluid: The weld pool has isotropic thermophysical properties, independent of strain rate.

Incompressible: The weld pool is assumed to be an incompressible liquid in that the dependence of its density onto both temperature and pressure is neglected. This assumption is argued for because any variations in the dynamic pressure and the temperature of the weld pool negligibly affect the density of the liquid metal. However, the thermal expansion of the fluid responsible for the buoyancy effect is assumed to behave according to the Boussinesq approximation. This means that a slight variation in density is accounted for in the gravitational effect, which allows for the capture of buoyancy effects [1].

Resistive magnetic fluid: The weld pool is assumed to be magnetically resistive in that the magnetic field is decoupled from the fluid flow and no magnetic advection occurs. This hypothesis is considered due to the low magnetic Reynolds number of the liquid metal Re m << 1.

An additional hypothesis is posed here and concerns the flow regime of the weld pool, and is detailed in the following:

Laminar flow: The weld pool is assumed to flow under a laminar regime and so a turbulence model is not implemented. This approach is valid under the assumption that the resolution of the flow, for the range of welding configurations considered in this thesis, is sufficiently fine. However, the claim that the flow remains laminar can be challenged because the Reynolds numbers in the simulated weld pool can reach up to Re ≳ 10 3 , which makes the flow regime difficult to define 5 . However, because a deeper investigation of the true flow regimes of different weld pools is out of the scope of this thesis, the laminar flow model is retained.

Having discussed the main hypotheses taken in the fluid model of the weld pool, the modelled mass and momentum conservation laws are presented in the following.

Conservation of mass:

The equation used to model the conservation of mass in the weld pool is as described below:

∇ • v = 0 (2.22)
where v ∈ Ω pool is defined as the velocity field in the weld pool. Interpreting equation (2.22), any variation in the volumetric flow in one direction is to be compensated by flow along the other two spatial dimensions.

Conservation of momentum:

The equation used to model the conservation of momentum in the incompressible weld pool is as follows:

ρv • (∇v) = -∇p + ∇•µ * ∇v + ∇ t v + f Lor + f Bou (2.23)
where p is the weld pool pressure, ρ is the pool density and is considered constant and µ * is the strengthened dynamic viscosity and is a function of temperature. Furthermore, the f Lor is the Lorentz as defined by equation (2.6), f Bou the buoyancy force. The strengthened dynamic viscosity µ * is used to strongly dampen the pool flow in the mushy zone. The forces are modelled in the following:

f Bou = -β ρ ref g(T -T ref ) (2.24)
where β is the thermal expansion coefficient, ρ ref and T ref the reference density and temperature values.

As for µ * , it is defined as the following:

µ * = µ f L + µ damp (1 -f L ) (2.25)
where µ damp is a dampening viscosity, chosen at µ damp = 250 kg•m -1 s -1 in this work and f L is the liquid fraction, calculated by:

f L =          0 T < T s T -T s T l -T s T s ≤ T ≤ T l 1 T > T l .
(2.26)

where the strengthened dynamic viscosity force is adapted from Voller [START_REF] Vaughan R Voller | An enthalpy method for convection/diffusion phase change[END_REF]. It serves to increase the dynamic viscosity of the weld pool by assigning to the pool the viscosity of the liquid metal µ * = µ when f L = 1 to µ * = A when f L = 0. The material properties are given in appendix A as they are reported for 316L steel by Pichler and an IAEA 6 report [START_REF] Pichler | Measurements of thermophysical properties of solid and liquid nist srm 316l stainless steel[END_REF][START_REF] Bobkov | Thermophysical properties of materials for nuclear engineering: a tutorial and collection of data[END_REF]. 5 As shown in chapter 4. 6 International Atomic Energy Agency

Boundary and interface conditions

With the mass and momentum conservation laws modelled as four partial differential equations, the definition of appropriate closure relations defined at the dim(Ω) -1 subdomains is necessary. Furthermore, because equations (2.20), (2.21), (2.22) and (2.23) determine both the kinematics and dynamics of the fluid subdomains, both kinematic and dynamic closure relations can be applied. The generalised kinematic and dynamic conditions on the arc-plasma and weld pool boundaries ∂Ω are presented first, following with the conditions at the shared fluid interface Γ between the arc-plasma and weld pool subdomains.

Boundary conditions for the arc-plasma subdomain

The arc-plasma subdomain Ω pla is bounded by an inlet flow kinematic condition at the nozzle, no-slip kinematic conditions at the rigid walls, zero stress dynamic boundary conditions at free boundaries and symmetry conditions 7 if symmetry axes or planes are defined. The generalised mass and momentum boundary conditions of the arc-plasma are defined as follows:

Inlet: u

∂Ω in = Vin ∥Ω in ∥ n (2.27)
Rigid wall: u

∂Ω wall = U wall (2.28)
Free boundary:

(n • σ • n) ∂Ω f ree = 0 (2.29)
Symmetry boundary:

     n • u ∂Ωsym = 0 (n • σ • t) ∂Ωsym = 0 (2.30)
where Vin is the inlet volume flow rate at the nozzle entry, n the normal vector, U wall the velocity of the rigid wall w.r.t to the laboratory frame of reference, σ the stress tensor and t = t 1 + t 2 is the total tangential vector. The stress tensor is defined as σ = -P I + µ ∇u + ∇ t u -2 3 ∇• uI .

Boundary conditions for the weld pool subdomain

The weld pool subdomain Ω pool is bounded by the solid boundary of the pool and a symmetry condition if symmetry axes or planes are defined. The generalised mass and momentum boundary conditions of the weld pool are defined as follows:

Solid boundary: v

∂Ω pool = U wall (2.31)
Symmetry boundary:

     n • v ∂Ωsym = 0 (n • σ • t) ∂Ωsym = 0 (2.32)
7 Generally consists of a mix of kinematic and dynamic boundary conditions.

Indeed, U wall is defined as the velocity of the rigid wall w.r.t the laboratory frame of reference; however, it is also used when modelling the displacement effects in a 3D welding configuration. Chapter 6 discusses the application in detail, where the displacement of the workpiece w.r.t the laboratory frame of reference (which is fixed on the cathode) allows for the capture of displacement effects for the steady state model used in this thesis. This modelling choice assigns the workpiece displacement velocity V displ to the wall boundaries that are in motion w.r.t the laboratory frame of reference U wall = V displ , and assigns a null velocity to the wall boundaries fixed at the laboratory frame of reference U wall = 0 m•s -1 .

Closure conditions at the arc-pool interface

Although both the interface Γ and boundary ∂Ω subdomains are defined at dim Ω -1 dimension, they require different closure conditions. Unlike boundary conditions, interface conditions require kinematic and dynamic relations that allow for a transmission of the mechanical quantities between the two subdomains or material phases. Applied to the model used in this thesis, the interface requires conditions that conserve mass when deformed by the superposed forces applied across it. However, the conservation laws, in their present form need reformulation when analysing interfaces. Indeed, equations (2.20), (2.21), (2.22) and

(2.23) are defined for an arbitrary element dΩ and thus requires mapping to an arbitrary interface element dΓ. The mapping is performed using the Gaussian pillbox method 8 . Note that in the following analyses, the ∇ s operator is the surface gradient operator. A formal derivation of the interface kinematics and dynamics is necessary because they manifest as the dominant driving mechanisms in the pool [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF][START_REF] Mills | Factors affecting variable weld penetration[END_REF][START_REF] Traidia | Multiphysics modelling and numerical simulation of GTA weld pools[END_REF].

Thus, at the fluid interface the conservation laws reduce to:

Kinematic conditions:

Starting from the general equation for mass conservation as stated in equation (2.19), it is integrated over an arbitrary Gaussian pillbox (visualised in 2D in figure 2.2) subdomain Ω that exists inside both the arc and pool subdomains Ω pla , Ω pool :

∂Ω (ρu) • n d∂Ω = δl 0 ((ρu) Γ -(ρu) Γ+dΓ ) • n Γ dl + Γ ((ρu) arc • n arc -(ρu) pool • n pool ) dΓ = 0 (2.34)
interested with the kinematics at the interface Γ, the thickness δl of the pillbox is assumed to be much smaller than the height H of the system, with δl << H. This transforms equation (2.34) into:

((ρu) Γ -(ρu) Γ+dΓ ) • n Γ • δl + Γ ((ρu) arc • n arc -(ρu) pool • n pool ) dΓ = 0 (2.35)
and reducing the thickness of the pillbox to zero, equation (2.35) becomes:

lim δl→0 ((ρu) Γ -(ρu) Γ+dΓ )•n Γ •δl+ Γ ((ρu) arc •n arc -(ρu) pool •n pool ) dΓ = 0 = Γ ((ρu) arc -(ρu) pool )•n dΓ (2.36)
reverting to the original notation and for an arbitrarily sized interface dΓ, equation (2.36) reduces to:

((ρu) arc -(ρu) pool ) • n → (ρ arc u -ρ pool v) Γ • n = 0 (2.37)
where the n is the normal vector at the interface Γ. The generalised mass conservation from equation (2.20) reduces to equation (2.37), and is aptly named the normal kinematic interface condition. Relating two unknown variables along an unknown direction, the condition requires a supplementary closure relation, which is derived from the level set equation [1] at the interface (the level set is visualised in figure 2.3), as shown in the following:

F = z -h z (x, y) = 0 (2.38)
where F is the level set, h z (x, y) the height of the deformed interface w.r.t the initial configuration Γ 0 (see figure 2.3). Furthermore, because the level set is null, its material derivative is also null, as described in the following:

dF dt = dz dt - dh z dt = 0 → dz dt = dh z dt (2.39)
transforming the material derivatives to their Eulerian descriptions:

dz dt = u • ẑ dh z dt = ∂h z ∂t + (u • x) ∂h z ∂x + (u • ŷ) ∂h z ∂y              substituting -------→ eq.(2.39) ∂F ∂t -u • ∇F = (u • ẑ) - ∂h z ∂t -u • ∇h z = 0 (2.40)
where the newly introduced unknown variable h z is related to the normal at the interface n by:

n = ∇F ∥∇F ∥ (2.41)
so substituting equation (2.41) into (2.40), the resulting equation becomes:

u • n = - ∂F ∂t 1 ∥∇F ∥ (2.42)
where for an interface at its steady state equilibrium position ∂F ∂t = 0 which transforms equation (2.42) into the impermeability condition:

u • n = 0 where u = (u or v) (2.43)
where the combination of equations (2.37) and (2.43) are the interface conditions for the normal velocity components at the interface. This leaves the tangential velocity components requiring closure, and so the no-slip hypothesis at fluid interfaces is introduced, which holds for most fluid systems modelled with continuum theory, and is assumed to hold here [START_REF] Nikolaos | Chapter 9 -boundary-layer flow[END_REF]:

u • t = v • t (2.44)
where additional closure equations defining the unknown velocity components equations (2.37), (2.43) and

(2.44) to the different forces at the interface are necessary.

Dynamic conditions:

Starting from the general equation for momentum conservation which resembles 9 that of equation (2.21), it is integrated over an arbitrary Gaussian pillbox (cf. figure 2.2) subdomain Ω that exists inside both the arc and pool subdomains Ω pla , Ω pool :

Ω ρu • ∇u dΩ = Ω volu. f Ω dΩ + Ω ∇ • σ dΩ + Γ inte. f Γ dΓ (2.45)
where the f Ω represents volume forces and f Γ surface forces at the interface. The different identified forces acting inside the pillbox control volume Ω are the following:

volu. f Ω = f g + f Lor inte.
f Γ = f γ 9 A term related to surface forces is added, because the control volume used in the integral formulation crosses the interface.

where f g is the gravitational force and f γ is the surface tension. To facilitate the analysis, each term in equation (2.45) is discussed individually first before the simplification and reformulation of the dynamic condition. First, the advection term is analysed for δl << H:

Ω ρu • ∇u dΩ divergence ------→ theorem ∂Ω ρu • (u • n) d∂Ω = Γ ρu • (u • n) arc -ρu • (u • n) pool dΓ + O(δl) (2.46)
where the O(δl) represents the integration terms of order δl. Next, the gravitational force term is reformulated as a geopotential force where the divergence theorem is then applied to it for δl << H, in:

f g = ρg = ρ ∇Φ where Φ = gz + cnst Ω ρ ∇Φ dΩ divergence ------→ theorem ∂Ω (ρΦ) • n d∂Ω = Γ (ρ arc -ρ pool )Φ Γ • n dΓ + O(δl) (2.47)
next, the Lorentz force term is analysed for δl << H:

Ω f Lor dΩ = Γ f Lor dΓ δl (2.48)
next, the stress flow term is analysed for δl << H:

Ω ∇ • σ dΩ divergence ------→ theorem ∂Ω σ • n d∂Ω = Γ (σ arc -σ pool ) • n dΓ + O(δl) (2.49)
next, the surface tension force term is analysed along the interface Γ as it is derived using either a variational or differential approach in [START_REF] Wu | Numerical analysis of both front-and back-side deformation of fully-penetrated gtaw weld pool surfaces[END_REF][START_REF] Tryggvason | Direct numerical simulations of gasliquid multiphase flows[END_REF]:

f γ = M arangoni f orce ∇ s γ + surf ace tension (γ∇ s • n)n (2.50)
where γ is the surface tension coefficient and is considered to only be a function of surface temperature and sulphur content in weld pools [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF][START_REF] Traidia | Multiphysics modelling and numerical simulation of GTA weld pools[END_REF][START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF]. The Marangoni force can be reformulated into its explicit temperature dependent form by following the mathematical procedure done between equations (2.16) and

(2.18), transforming it into:

∇ s γ = ∂γ ∂T ∇ s T (2.51)
where γ as a function of temperature and the concentration of sulphur c S in parts per million [ppm] is described by Sahoo's law [START_REF] Sahoo | Surface tension of binary metal surface active solute systems under conditions relevant to welding metallurgy[END_REF]:

γ(T, c S ) = γ f -A g (T -T f ) -RT Γ s ln 1 + k 1 c S exp - ∆H 0 RT (2.52) ∂γ ∂T (T, c S ) = -A g -RΓ s ln 1 + k 1 c S exp - ∆H 0 RT - ∆H 0 exp - ∆H 0 RT c S Γ s 1 + exp - ∆H 0 RT c S T (2.53)
where γ f is the surface tension of pure iron at its temperature of fusion T f , A g = -∂γ ∂T (T, 0), T the pool surface temperature, R the ideal gas constant, Γ s the surface excess at saturation, k 1 a segregation entropy constant, ∆H 0 the enthalpy of adsorption. The operations performed above, reassembled together become:

Γ ρu • (u • n) arc -ρu • (u • n) pool dΓ = Γ (ρ arc -ρ pool )Φ • n dΓ + Γ (σ arc -σ pool ) • n dΓ + Γ f γ dΓ + Γ f Lor dΓ δl + O(δl)
and once again, the thickness of the Gaussian pillbox is taken to zero by applying the lim δl→0 to the above equation and taking an arbitrarily sized interface dΓ, the pillbox momemtum equation (2.45) maps to the following dynamic interface condition:

ρu • (u • n) arc -ρu • (u • n) pool = (ρ arc -ρ pool )gh z • n + (σ arc -σ pool ) • n + ∂γ ∂T ∇ s T + (γ∇ s • n)n (2.54)
where both the normal and tangential dynamics are included. It is crucial to remark that the advection terms at the left hand side of equation (2.54) invoke the kinematic interface relation from equation (2.37).

The advection terms in the dynamic interface condition are thus critical if there are sources or sinks in the mass balance of the system; where in TIG welding the sources and sinks can be related to both phase change and/or a dumping of filler material into the pool subdomain. However, since this thesis is limited to studying a TIG system without evaporation effects nor a filler metal, phase change and filler dumping cases are out of the scope of this work and their analysis is reserved for a future study. Thus, using the impermeability condition in equation (2.43), this simplifies equation (2.54) to:

0 = (ρ arc -ρ pool )gh z • n + (σ arc -σ pool ) • n + ∂γ ∂T ∇ s T + (γ∇ s • n)n (2.55)
which will be used as the dynamic boundary condition in this thesis. Both the arc shear stress and the Marangoni effect at the interface are dominant driving forces for the pool hydrodynamics [START_REF] Mills | Factors affecting variable weld penetration[END_REF].

Interface hypotheses

A brief description of the two interface hypotheses treated in this thesis are presented. The fixed interface hypothesis assumes that the API remains fixed while the deformable interface hypothesis assumes that the API is free to deform.

Fixed interface:

The fixed interface hypothesis maintains the arc-pool interface at its initial configuration Γ 0 by assuming that the sum of normal reaction forces at dΓ 0 is at equilibrium. The sum of forces at fixed interface element is schematised in figure 2.4, where neither the surface tension force nor the gravity force act on dΓ. This is because h z = 0 which makes null the gh z and γ∇ s • n terms. Thus for fixed interfaces, the normal arc and pool stress forces are at equilibrium while the sum of the tangential arc and pool stress forces are equal to the Marangoni force. Summarising the kinematics and dynamics at dΓ, the set of interface conditions to interface conditions to be solved at the deformable interface are:

kinematic impermeability condition: (u -v) • n = 0 ; v • n = 0 kinematic no-slip condition: (u -v) • t = 0 level set: v • ẑ = v • ∇h z ; n = ∇F ∥∇F ∥ ; F = z -h z dynamic condition: 0 = (ρ arc -ρ pool )gh z • n + (σ arc -σ pool ) • n + ∂γ ∂T ∇ s T + (γ∇ s • n)n

Energy conservation law

The TIG welding system is made up of three bodies that transform electrical energy into heat and that transfer the generated heat to the workpiece, which melts to form the liquid weld pool. The cathode, arc-plasma and the workpiece subdomains are simultaneously subject to electro-thermal, heat transport and phase change phenomena. Thus, to effectively capture the range of thermal effects in the TIG system, the energy conservation law must be modelled in a manner compatible with both conjugate heat transfer and phase change phenomena. To this end, energy conservation in both the cathode and arc-plasma are modelled using the temperature variable as the primary unknown variable, while energy conservation in the workpiece is modelled using the enthalpy variable. The mixed variable method is used because the temperature based models allow for a simpler expression of the electro-thermal phenomena in the cathode and arc-plasma, while the enthalpy based model is adapted for the capture of phase change effects in the workpiece. The models are briefly presented in the following and are based on the theses of Brochard and

Nguyen [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF].

Temperature based heat transfer model

The cathode and arc-plasma interact electro-thermally to transmit and transform electrical energy into heat energy useful for TIG welding. To capture these effects, temperature based models in both the cathode and arc-plasma subdomains are set up. The cathode is modelled as a solid because for most TIG welding time scales, the cathode remains rigid and does not undergo fusion [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Mills | Factors affecting variable weld penetration[END_REF][START_REF] Tanaka | An introduction to physical phenomena in arc welding processes[END_REF][START_REF] Chihoski | The rationing of power between the gas tungsten arc and electrode[END_REF]. As for the arc-plasma, it is hypothesised to maintain LTE even though its colder edges (particularly the electrode boundaries) deviate from LTE [START_REF] Gleizes | Thermal plasma modelling[END_REF][START_REF] Lago | A numerical modelling of an electric arc and its interaction with the anode: Part I. the two-dimensional model[END_REF][START_REF] Lowke | LTE diffusion approximation for arc calculations[END_REF]. To simplify modelling of the deviation from LTE, equivalent electro-thermal source and sink terms are added at both the CPI and API. Thus, the arc-plasma is modelled as a fluid governed by the LTE hypothesis. The modelled equivalent sources and sinks at the electrode-arc interfaces are adapted from Brochard's thesis [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF].

Cathode subdomain

The energy conservation law in the cathode subdomain is from Brochard [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF], and is modelled by the following equation:

0 = ∇• λ∇T + s Joule (2.56)
where the thermal conductivity λ as a function of temperature for thoriated tungsten is given in appendix A. Heat in the cathode is generated by the Joule effect as modelled by equation (2.5).

Arc-plasma subdomain

The energy conservation law in the arc-plasma is from Brochard [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF], and is modelled by the following equation:

ρc p u•∇T = ∇• λ∇T + s Joule -s Ray,pla (2.57)
where the mass density ρ, the specific heat c p and the thermal conductivity λ of argon are given in appendix A. Heat generation in the arc-plasma is modelled by the Joule effect as it is given by equation (2.5), while the total radiated heat is modelled by the emissive volumic radiation term s Ray,pla as adapted from Lago et al. [START_REF] Lago | A numerical modelling of an electric arc and its interaction with the anode: Part I. the two-dimensional model[END_REF]. The volume radiation term is modelled in the following manner:

s Ray,pla = 4πϵ n (2.58)
where ϵ n is the total argon plasma emissivity term and is given as a function of temperature in appendix A, as it was reported by Lago et al. [START_REF] Lago | A numerical modelling of an electric arc and its interaction with the anode: Part I. the two-dimensional model[END_REF].

Enthalpy based heat transfer model

The heat transported by the arc-plasma melts the weld zone of the workpiece, inducing phase change over the course of a welding operation. The phase change process is generally strongly temperature dependent and nonlinear in that the enthalpy of the studied material rises significantly over a small temperature range 10 . Thermodynamically the enthalpy h is related to the temperature state variable T and the pressure state variable p by the following relation [START_REF] Kurz | Fundamentals of Solidification: Fourth Revised Edition[END_REF]:

h = T T ref c p dT + p p ref ∂h ∂p dp (2.59)
where T ref and p ref are some arbitrary reference temperature and pressure. The dependence of the enthalpy onto the pressure is negligible for the weld pool system because any dynamic variations in the pool pressure are asymptotically small w.r.t its thermodynamic pressure [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF][START_REF] Kurz | Fundamentals of Solidification: Fourth Revised Edition[END_REF]. This allows for a simplified isobaric enthalpy relation, which transforms equation (2.60) into:

if h = f (T ) → ∂h ∂p ≈ 0 ; h = T T ref c p dT (2.60)
which applies for both the solid and liquid states of the weld material. The notion of enthalpy is important when discussing phase change problems because it masks the strongly nonlinear dependence of phase change onto temperature by acting as an integral quantity. This can be understood by comparing the expression of the specific heat as a function of temperature to the expression of the enthalpy as a function 10 Or, in the case of pure metals, the enthalpy rises over an isotherm. See appendix B for a detailed description.

of temperature, where the specific heat for a material undergoing fusion is:

c p (T ) =            c ps for T < T s L f 1 2π(∆T m /a) 2 exp(- (T -T c ) 2 2(∆T m /a) 2 ) for T s ≲ T ≲ T l c p l for T < T l (2.61)
where T s and T l are the solidus and liquidus temperatures of the material, L f the latent heat of fusion,

T m = T l -T s the mushy temperature range, T c = T l +Ts 2
the mushy center temperature and a a sharpness factor to be > 1. For details on how the specific heat is modelled, see appendix B. While the enthalpy for a material undergoing fusion is:

h(T ) =        T T ref c ps dT for T < T s Ts T ref c ps dT + L f f L for T s ≤ T ≤ T l Ts T ref c ps dT + T T l c p l dT + L f for T > T l (2.62)
where f L is defined as the liquid fraction from equation (2.26). Indeed, in the range of the mushy zone where phase change occurs T s ≤ T ≤ T l , the model used for the specific heat invokes temperature gradients that are significantly stronger than those of the enthalpy model. This is because the enthalpy variable intrinsically contains the latent heat contributions in its formulation, thus masking the small temperature scales of the specific heat variable which introduce strong non-linearities to temperature based heat equations. The differences in the temperature scales of equations (2.61) and (2.62) are schematised in figure 2.6.

Thus, using the enthalpy variable instead of the temperature variable in the energy conservation model is

ρv•∇h = ∇• λ∇ h c p + s Joule (2.63)
where ρ, v, c p and λ are the density, velocity, specific heat and thermal conductivity of the weld pool, respectively. The material properties of the studied workpieces are presented in appendix A as they are reported by Brochard and Pichler respectively [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Pichler | Measurements of thermophysical properties of solid and liquid nist srm 316l stainless steel[END_REF]. Note that v is equal to the workpiece displacement velocity V displ w.r.t the laboratory reference frame in the solid workpiece subdomain. Moreover, the heat source by the Joule effect is modelled by equation (2.5). Equation (2.63) is based on the model used in Nguyen's thesis [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF].

Boundary and interface conditions

With the energy conservation laws modelled as three partial differential equations, the definition of appropriate closure relations defined at the dim(Ω) -1 subdomains is necessary. Furthermore, because equations (2.56), (2.57) and (2.63) are defined in the volume subdomains of the system appropriate temperature continuity and flux interface conditions are needed to account for the heat transfer phenomena across the interfaces. The generalised thermal conditions on the cathode, arc-plasma and workpiece boundaries ∂Ω are presented first, following with the condtions at the shared interfaces Γ between the cathode and arc-plasma, and the arc-plasma and workpiece subdomains.

Boundary conditions for the cathode and arc-plasma subdomains

With both the cathode and the arc plasma modelled with the temperature variable, their boundary conditions are similarly applied. The cathode and arc-plasma are bounded by fixed temperature, free boundary flux (also a thermal symmetry condition) and radiation conditions, which are applied using the Stefan-Boltzmann radiation model [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Theodore L Bergman | Fundamentals of heat and mass transfer[END_REF]. The generalised thermal boundary conditions of the cathode and arc-plasma are defined as follows:

Fixed temperature:

T ∂Ω f ixed = T f ixed (2.64)
Free boundary or symmetry:

(-λ∇T • n) ∂Ω f ree = 0 (2.65) Radiation: (-λ∇T • n) ∂Ω rad = ϵσ B (T 4 -T 4 ∞ ) (2.66)
where T f ixed is a chosen temperature that remains fixed throughout a simulation, ϵ, σ B and T ∞ the emissivity, Stefan-Boltzmann constant and the free stream temperature, respectively, and are given in appendix A.

Boundary conditions for the anode subdomain

With the workpiece (or anode) modelled with the enthalpy variable, its temperature dependent boundary conditions are specified using temperatures that were mapped from their associated enthalpy variable (see section 3.5.1 for details). The workpiece is bounded by fixed enthalpy, free boundary flux (also a thermal symmetry condition), radiation and convection conditions. The radiation is similarly modelled using the Stefan-Boltzmann model while Newton's law of cooling is used to model the convection condition [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Theodore L Bergman | Fundamentals of heat and mass transfer[END_REF].

The generalised thermal boundary conditions of the anode are defined as follows:

Fixed temperature: h

∂Ω f ixed = h(T f ixed ) (2.67)
Free boundary or symmetry:

(-λ∇T • n) ∂Ω f ree = 0 (2.68) Radiation: (-λ∇T • n) ∂Ω rad = ϵσ B (T 4 -T 4 ∞ ) (2.69) Convection: (-λ∇T • n) ∂Ωconv = h conv (T -T ∞ ) (2.70)
where h conv is the convection coefficient.

Closure conditions at the interfaces

Unlike boundary conditions, interface conditions require both temperature continuity and flux relations that allow for a transmission of the thermal quantities between the two subdomains or material phases.

The thermo-electric sources and sinks used to model the electrode boundaries at the interfaces require conditions that conserve energy when superposed fluxes are applied across it. However, the conservation laws, in their present form need reformulation when analysing interfaces. Indeed, equations (2.56), (2.57),

(2.63) are defined for an arbitrary volume element dΩ and thus require mapping to an arbitrary interface element dΓ. The closure conditions at the CPI and API are presented in the following.

Cathode-plasma interface:

First, temperature continuity at the interface is assumed. The continuity of the thermal field is valid at the interface of perfectly conjugate materials [START_REF] Theodore L Bergman | Fundamentals of heat and mass transfer[END_REF]; which, based on the LTE hypothesis used in this work, is the case for the CPI. This is mathematically expressed as:

T Γcat = T Γarc (2.71)
Next, the flux condition is expressed by applying the Gaussian pillbox method as was previously demonstrated, which transforms the energy conservation law (equation (2.57)) at the CPI to the following form:

0 = (q arc -q cat ) • n + inte. s Γ CP I (2.72)
where q = -λ∇ s T and s Γ represents the heat sources/sinks at the interface. Note that the advection terms reduce to null because the CPI is a rigid undeformable interface. The identified sources, as modelled by Brochard [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF], are:

inte.

s Γ CP I = s Neu -s Emi -s SB
where s Neu is the ionic neutralisation source term, s Emi and s Ray,cat are the thermoionic emission and Stefan-Boltzmann (see equations (2.66), (2.69)) sink terms. The electric heat terms are functions of electric current density, as presented in the following:

s Neu = j Ion V i s Emi = j e -W cat
where V i and W cat are the ionisation potential and the cathode work function, respectively. The ionic current j Ion and emitted current j e -are calculated in the following manner:

j e -= j CP I -j Ion j RD = A r T 2 exp -W ef f k B T j Ion = max(j CP I -j RD , 0) j CP I = σ * ∇ϕ • n
where A r is the Richardson constant, W ef f the effective work function of the cathode and k B Boltzmann's constant and they are given in appendix A. Thus, equation (2.72) becomes:

0 = (q arc -q cat ) • n + j Ion V i -j e -W cat -ϵσ B (T -T ∞ ) 4 (2.73)

Anode-plasma interface:

Similarly to the CPI, the continuity of the thermal field based on the LTE hypothesis is first addressed, where temperature continuity, using both the temperature and enthalpy variables becomes:

T Γarc = h h ref 1 c p (T ) dh Γano + T ref (2.74)
where h ref is an arbitrary reference enthalpy calculated at some reference temperature T ref .

Next, the flux condition is expressed by applying the Gaussian pillbox method as was previously demonstrated, which transforms equation (2.63) into:

(ρu • c p T ) • n arc -(ρv • h) • n ano = (q arc -q ano ) • n + inte. s Γ AP I (2.75)
It is once more crucial to remark that the advection terms at the left hand side of equation (2.75) invoke the kinematic interface relation from equation (2.37). The advection terms in the flux condition are thus critical if there are sources or sinks in the mass balance of the system. With no mass sources nor sinks considered in this work, the left hand side is null. This simplifies equation (2.75) to:

0 = (q arc -q ano ) • n + inte. s Γ AP I (2.76)
where the identified source terms are:

inte.

s Γ AP I = s Abs -s SB
where s Abs is the electronic absorption term and s SB the Stefan-Boltzmann radiation term (see equations

(2.66), (2.69)). The electronic absorption term is calculated in the following manner:

s Abs = j AP I W ano for j AP I = -σ * ∇ϕ • n
where W ano is the work function of the anode and is given in appendix A. Thus equation (2.76) becomes:

0 = (q arc -q ano ) • n + j AP I W ano -ϵσ B (T -T ∞ ) 4 (2.77)
where q arc = -λ∇T and q ano = -λ∇ h cp , thus allowing for an interface relation consistent with the mixed variable approach of this section.

Summary

The modelled electrodynamics and conservation laws that capture the dominant phenomena at play in the TIG system were presented in this chapter. The models were discussed as they apply for volume elements, while their mappings to the cathode-arc and anode-arc interfaces were derived. Furthermore, the relevant generic boundary conditions for each model were detailed.

The used electrodynamic model respects the simplified LTE and continuum hypotheses. This model captures the dominant electromagnetic phenomena that occur at the mesoscopic scale. This approach allows, in a simplified manner, to calculate both the Joule effect and the Lorentz force in the TIG system. This is crucial because both the Joule effect and the Lorentz force are the dominant heat and momentum sources in the arc-plasma.

The models used for the momentum and mass conservation laws in both the arc-plasma and the weld pool allow for the resolution of the dominant fluidic effects which are important for the capture of thermal effects in the global model. Both the kinematics and dynamics at the arc-pool interface are formally derived and presented, because they manifest as the dominant driving mechanisms in the pool. Therefore, ensuring that both the arc shear and Marangoni effect are properly manifest in the dynamic interface condition is paramount. Furthermore, the generalised interface conditions bring into context the importance of nonnegligible advective terms that can govern the interface thermo-hydraulics if source/sink terms are added to the system.

The modelling of the energy conservation law in both the arc-plasma and the workpiece was presented in a manner compatible with both conjugate heat transfer and phase change phenomena. Effectively, a mixed variable formulation was presented in which the temperature variable was used for the cathode and arc-plasma model and the enthalpy variable for the workpiece model. Furthermore, the interface conditions that incorporate the equivalent electronic heat sources at the interfaces in a manner consistent with the mixed variable method were presented.

In this chapter the relevant physical phenomena pertaining to the TIG process were presented in their mathematical form. To solve the equations of the global model, numerical methods, coupling techniques and algorithms are needed. The development and implementation of these methods are detailed in the following chapter.

Chapter 3

Numerical methods, coupling techniques and algorithms

The beauty of the universe is witnessed through numbers...

Geometric preliminaries

To facilitate the numerical discussion performed in this chapter, key geometric aspects are presented and are based on the geometric discussion of section 2.2 of the previous chapter. The Γ AP I , being a subdomain where two fluids interface, is further decomposed into two subdomains. Thus, Γ AP I is numerically treated in the following manner:

Γ AP I = Γ plaapi ∪ Γ anoapi
where Γ plaapi and Γ anoapi are geometrically overlapping subdomains that are numerically stored as separate objects in the computer program. These different numerical objects make up the interfacial nodes of the mesh of their respective subdomains Ω pla and Ω ano . A simplified representation of the geometric structure of the system is presented in figure 3.1. The meshing method of overlapping the interfacial nodes at Γ AP I is also represented where the geometrically overlapping nodes can be assigned to different arbitrary functions.

The schematised functions

f = [f 1 , f 2 , f 3 , ..., f N ] and g = [g 1 , g 2 , g 3 , .
.., g N ] represent any mix of functions defined at the nodes at either side of the interface, that can have some linear relationship. The linear constraint can be prescribed using the Lagrange Multiplier technique [START_REF] Pellet | Dualisation des conditions aux limites[END_REF]. Moreover, the interfacial nodes at either the Γ plaapi or the Γ anoapi side can be loaded with surfacic source or sink terms.

A brief discussion of the chosen numerical methods

This section presents the basics of the chosen computational methods, used to solve for different TIG welding configurations. A brief discussion of the Galerkin FEM method used to solve for the electric, thermal, momentum and mass equations is presented, followed by the Least Squares FEM (LSFEM) method used to solve the magnetic equations. Finally, the different stabilisation methods used for certain mathematical operators are briefly presented.

The Galerkin FEM approach

In an attempt to be as concise as possible, the details of the Galerkin FEM discretisation and the variational weighting of the residual vector are assumed trivial. This allows for a more direct discussion of the way the residual vectors are constructed for each of the partial differential equations (PDE) of the physical model. The interested reader is referred to the texts by Reddy [START_REF] Reddy | An Introduction to the Finite Element Method[END_REF] and Kuzmin [START_REF] Kuzmin | Finite element methods for computational fluid dynamics: a practical guide[END_REF] for the theoretical details of the Galerkin FEM method. Thus, we can consider that the general form of the weak variational formulation of the residual vector R ψ is as it appears in the following algebraic system for some variable ψ, defined over an arbitrary mesh:

{R ψ } = 0 = [M ]{ψ} -{F } (3.1)
where [M ] and {F } are the total matrix and the force vector, respectively. The total matrix is the discrete, variational form of the differential operators and can be made up of the advection, convection, gradient and/or Laplacian operators. Depending on the discretised PDE, the total matrix can be a function of the primal variable ψ itself, and/or can be dependent on temperature. In fact, the modelled PDEs in this thesis are strongly non-linear, and so a linearisation step can be necessary before solving for ψ; this is discussed in section 3. The electric model, based on equation (2.1), in its discrete variational form, manifests as:

[M ]{ϕ} -{F } = 0 = {R} ϕ ≡ Ω (σ * ∇ϕ) • (∇N ϕ ) dΩ - ∂Ω imp (j imp • n)(N ϕ ) d∂Ω imp (3.2)
where (a) in Ω cat :

the primal variable is the temperature variable, {T }, and the residual looks like:

0 = R T = Ωc (λ∇T ) • (∇N T ) dΩ c - Ωc (σ∇ϕ) • ∇ϕ (N T ) dΩ c - Γ cpi j Ion V i -j e -W cat -ϵσ B (T 4 -T 4 ∞ ) N T dΓ cpi (3.4) (b) in Ω pla :
the primal variable is the temperature variable, {T }, and the residual looks like:

0 = R T = Ωp (ρc p u∇T )N T dΩ p + Ωp (λ∇T ) • (∇N T ) dΩ p - Ωp (σ∇ϕ) • ∇ϕ (N T ) dΩ p + Ωp (4πϵ n )N T dΩ p (3.5) (c) in Ω ano :
the primal variable is the enthalpy variable, {h}, and the residual looks like:

0 =R h = Ωa (ρv∇h)N h dΩ a + Ωa (λ∇T ) • (∇N h ) dΩ a - Ωa (σ∇ϕ) • ∇ϕ (N h ) dΩ a - Γ api (j AP I W ano ) -ϵσ B (T 4 -T 4 ∞ ) N h dΓ api - ∂Ωconv h conv (T -T ∞ )N h d∂Ω conv (3.6)
To stabilise the advection terms in the numerical model, a Streamline Upwind Petrov-Galerkin (SUPG) [START_REF] Thomas | Recent progress in the development and understanding of supg methods with special reference to the compressible euler and navier-stokes equations[END_REF] numerical diffusion term is added to the residuals. However, for brevity the term is not explicitly shown. Furthermore, the linearisation method and the treatment of the conjugate heat transfer problem with different primals {T } and {h} is discussed in sections 3.4 and 3.5.1, respectively.

Momentum and mass transfer model:

The momentum and mass transfer model, based on equations (2.20), (2.21), (2.22) and (2.23), modifies system (3.1), to manifest as:

[M (ρ, u, v, µ)]{u, P, v, p} = {F } (3.7)
where system (3.7) combines the momentum and mass conservation equations and is domain dependent, but is undefined for the Ω cat domain as it is modelled as a solid. Thus, nonlinear system (3.7) depends on Ω pla and Ω ano in the following manner:

(a) in Ω pla :

the primal variables are the velocity and pressure variables, {u} and {P }, and the residual looks like:

0 = R u = Ωp (ρu • ∇u) • (N u ) dΩ p + Ωp µ ∇u + ∇u t : ∇N u + ∇N t u dΩ p - Ωp 2 3 µ(∇ • u)(∇ • N u ) dΩ p - Ωp (-σ∇ϕ) ∧ B • N u dΩ p - Ωp P (∇ • N u ) dΩ p (3.8) 0 = R P = Ωp (∇ • u)N P dΩ p + Ωp d(ln ρ) dT ∇T • u N P dΩ p (3.9) (b) in Ω ano :
the primal variables are the velocity and pressure variables, {v} and {p}, and the residual looks like:

0 = R v = Ωa (ρv • ∇v) • (N v ) dΩ a + Ωa µ ∇v + ∇v t : ∇N v + ∇N t v dΩ a - Ωa (-σ∇ϕ) ∧ B • (N v ) dΩ a - Ωa ρgβ(T -T ref ) • N u dΩ a - Γ AP I ( dγ dT ∇ s T ) • (N v ) dΩ a - Ωa p(∇ • N v ) dΩ a (3.10) 0 = R p = Ωa (∇ • v)N p dΩ a (3.11)
To construct a stable velocity-pressure coupling scheme in the above residuals, either: i) the chosen vector spaces for the velocity and pressure variables respect the Ladyzhenskaya-Babuska-Brezzi condition [START_REF] Kuzmin | Finite element methods for computational fluid dynamics: a practical guide[END_REF][START_REF] Reddy | An Introduction to the Finite Element Method[END_REF], where the pressure variable has a discretisation order smaller than that of the velocity variable, where for any mesh element:

u h = m i u i N u,i P h = n i P i N P,i (3.12)
where N u,i and N P,i are polynomial functions of different order where n < m; n + 1 = m;

ii) or they are stabilised with a polynomial pressure projection term based on the work or Dohrmann [START_REF] Clark R Dohrmann | A stabilized finite element method for the stokes problem based on polynomial pressure projections[END_REF].

These velocity-pressure considerations are necessary for the stability of the Stokes problem in the U -P formulation [START_REF] Kuzmin | Finite element methods for computational fluid dynamics: a practical guide[END_REF]. Furthermore, the advection terms are stabilised with a SUPG term [START_REF] Thomas | Recent progress in the development and understanding of supg methods with special reference to the compressible euler and navier-stokes equations[END_REF].

For brevity, the stabilisers are not explicitly shown in the residuals. Furthermore, the linearisation method and the treatment of the fluid coupling method between the arc-plasma and the weld pool is discussed in sections 3.4 and 3.6, respectively.

The LSFEM approach

The LSFEM method is a technique generally adapted to the resolution of 1 st order PDEs, unlike Galerkin methods. Effectively, using the Galerkin FEM to solve 1 st order PDEs, without adding numerical diffusion can generate non-symmetric matrices with odd-even decoupling, and this leads to oscillatory solutions.

Thus, this motivates the use of the LSFEM technique to solve the magnetostatic model used in this thesis [START_REF] Reddy | An Introduction to the Finite Element Method[END_REF][START_REF] Jiang | The least-squares finite element method: theory and applications in computational fluid dynamics and electromagnetics[END_REF]. Motivating the use of LSFEM, a brief discussion of its advantages are first discussed and then the basics of the method are presented.

Advantages of LSFEM for the magnetostatic model

The magnetostatic model used in this thesis, benefits from the LSFEM method for a number of reasons, and the most relevant are discussed below:

• Gauss' Law of magnetism, modelled with the ∇ • B = 0 equation is easily incorporated into the LSFEM formulation of the magnetostatic model. In fact, assuming the redundancy of Gauss' law 3and thus its omission from a magnetic model, can generate parasitic stationary solutions. Thus, the ease of incoporation of Gauss' law in the LSFEM formulation makes the technique attractive [START_REF] Jiang | The least-squares finite element method: theory and applications in computational fluid dynamics and electromagnetics[END_REF].

• Although magnetostatic models can use the vector potential A and thus be solved for with the Galerkin FEM approach (see Nguyen's thesis for details [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF]); it can be suboptimal if the system is strongly sensitive to the magnetic field B. This is due to B = ∇ × A , thus any approximated solution:

A h ∝ A + O(δr m ) yields: B h ∝ ∇ × A + O(δr m-1 )
after differentiation using the ∇× curl operator. O(δr m ) is the discretisation error of order δr, which is the characteristic element size in a mesh. The reduction in the order of precision, going from m → m -1, can influence the precision of magnetic field dependent calculations [START_REF] Jiang | The least-squares finite element method: theory and applications in computational fluid dynamics and electromagnetics[END_REF]. For example, the velocity fields of TIG arcs are very sensitive to the Lorentz force, and thus can be consequently very sensitive to the precision of the magnetic field.

• The use of LSFEM to solve the magnetostatic model generates positive definite and symmetric matrices, which can be easily treated by linear solvers. Moreover, the generated linear systems, being stable and non-oscillatory do not require artificial diffusion methods for stabilisation, as is necessary for advection terms (see section 3.3.1) [START_REF] Reddy | An Introduction to the Finite Element Method[END_REF][START_REF] Jiang | The least-squares finite element method: theory and applications in computational fluid dynamics and electromagnetics[END_REF].

Basics of the method

Here, the basics of the LSFEM technique are introduced in as concise a manner as possible. The interested reader is referred to the text by Jiang [START_REF] Jiang | The least-squares finite element method: theory and applications in computational fluid dynamics and electromagnetics[END_REF] for more details. The technique consists of treating 1 st order boundary value problems, in the following manner:

M u = f (3.13)
where M is a 1 st order differential operator:

M u = n d n=1 M n ∂u ∂x n + M 0 u
looking for approximated solutions, and assuming that f ∈ L 2 (Ω) so that M : U → L 2 (Ω) and where the primal variable u is approximated by u h ∈ U :

u h = M i u i N i (3.14)
where the trial function N i ∈ U , and M is the number of nodes in an arbitrary element. Applying equation (3.14) to equation (3.13), the approximated equation becomes:

R = M u h -f 0 unless u h = u, the residual R 0.
The LSFEM method, unlike the Galerkin method, looks not to nullify the weighted residual Ω R i N j dΩ = 0, rather looks to minimize the squared distance between M u h and f , in:

I(u h ) = ∥R∥ 2 0 = Ω (M u h -f ) 2 dΩ ⩾ 0 (3.15)
where I(u h ) is defined as the functional. A necessary condition for u h to be a root imposing the minimum of the functional in equation 3.15 is that any perturbation to I(u h ) vanishes at u h , in:

lim δv→0 ∂ ∂δv I(u h + δvN i ) = 2 Ω (M N i ) T (M u h -f ) dΩ = 0
The minimization of the I(u h ) functional thus gives an algebraic system of the form:

Ω (M N j )(M N i )u i dΩ = Ω (M N j )f dΩ (3.16)
Equation (3.16) is the general form that results from an LSFEM discretization. Thus, the residual to the magnetostatic equations (2.2) and (2.3) can be expressed as:

0 = R B = Ω (∇ ∧ B) • (∇ ∧ N B ) dΩ + Ω (∇ • B)(∇ • N B ) dΩ - Ω (-σ∇ϕ) • (∇ ∧ N B ) dΩ (3.17)

The multiphysics computational toolbox

The so called multiphysics computational toolbox, as it is defined in this work, is a small kit of coupling

Coupling techniques

The coupling techniques discussed in this subsection are by no means comprehensive, and serve only to briefly present the methods used in this thesis. The interested reader is referred to Patricot's [START_REF] Patricot | Couplages multi-physiques : évaluation des impacts méthodologiques lors de simulations de couplages neutronique/thermique/mécanique[END_REF] and Viot's [START_REF] Viot | Couplage et synchronisation de modeles dans un code scenario d'accidents graves dans les reacteurs nucleaires[END_REF] theses, which offer more generalised presentations of the mathematical techniques used to couple different modules together. The coupling techniques used in this work are either a variant, and/or a combination of the Gauss-Seidel and the Newton-Raphson methods. The two methods are presented in the following:

1. Gauss-Seidel:

The Gauss-Seidel method is a simple method, easy to implement, and is a method that weakly cou-ples 5 two (or more) modules together. This coupling technique is widely used in the literature, where once the coupling is ensured and the mathematical variables are properly transmitted, the coupled system is iterated until convergence is achieved 6 .

With k representing the iteration number; and starting from an initial guess ψ 0 1 , ψ 0 2 , . . . , ψ 0 l the equations are iterated in the following manner, for l different modules:

             ψ k 1 = f 1 (ψ k-1 1 , ψ k-1 2 , . . . , ψ k-1 l ) ψ k 2 = f 2 (ψ k 1 , ψ k-1 2 , . . . , ψ k-1 l ) . . . ψ k l = f l (ψ k 1 , ψ k 2 , . . . , ψ k-1 l ) (3.18)
The stopping criterion is generally based on increment size where

ϵ k 1 = ψ k 1 -ψ k-1 1 , ϵ k 2 = ψ k 2 -ψ k-1 2 , . . . , ϵ k l = ψ k l -ψ k-1 l
, and convergence is assumed 7 once ϵ k l → 0. Coupling algorithm (3.18) generally converges to a solution if the initial guesses ψ 0 1 . . . ψ 0 l are well chosen and within the vicinity of the roots of the algebraic system [START_REF] Patricot | Couplages multi-physiques : évaluation des impacts méthodologiques lors de simulations de couplages neutronique/thermique/mécanique[END_REF].

Newton-Raphson:

The Newton-Raphson method is a classic method which is commonly used to solve non-linear and coupled systems. The method has many variants based on approximating the original algorithm, and can be generally classified as Quasi-Newton methods [START_REF] Viot | Couplage et synchronisation de modeles dans un code scenario d'accidents graves dans les reacteurs nucleaires[END_REF]. The approach is based on a linearisation of the coupled system of equations and iterating over the generated linear systems until convergence is achieved. If no approximations to the method are undertaken, the coupling is classified as globally monolithic or strong, where the algebraic systems and their couplings are solved simultaneously. The original method applied to a multiphysics model is applied to an algebraic system of equations similar to equation (3.1), in the following manner:

for any global set of variables ψ m , operators M mn (ψ n ) and duals F m , for l different modules,

   F 1 . . . F l    -    M 11 . . . M 1l . . . . . . . . . M l1 . . . M ll       ψ 1 . . . ψ l    =    R 1 . . . R l    =    0 . . . 0    (3.19)
and solving for the roots of system (3. [START_REF] Cuvelier | Some numerical methods for the computation of capillary free boundaries governed by the navier-stokes equations[END_REF]) is in most cases impossible without linearising it. Applying a 1 st order Taylor Expansion to the multivariable system around the roots ψ sol of system (3.19):

{R(ψ sol + ψ k )} = 0 ≈ {R(ψ k-1 )} + [J ψ ] k-1 {ψ sol -ψ k }
applying recursively the 1 st order Taylor Expansion a k number of times, assuming an initial guess of {ψ 0 } ≈ {ψ sol } -{δψ k } -{δψ k-1 } . . . is not too far away from the roots {ψ sol } of system (3.19), 5 Weak coupling refers to a coupling where the equations are not coupled and solved simultaneously in the same algebraic system. 6 Or stopped based on a mathematical criterion. 7 For badly conditioned systems of equations, the residual of each modelled equation is also required to tend to 0.

the linearised system, at iteration k, to solve becomes:

[J ψ ] k-1 {δψ} k = -{R ψ } k-1 (3.20)
where {δψ k } = {ψ k } -{ψ k-1 }, and where the global Jacobian matrix [J ψ ] is defined as:

[J ψ ] =    T 11 . . . T 1l . . . . . . . . . T l1 . . . T ll    =        ∂R 1 ∂ψ 1 . . . ∂R 1 ∂ψ l . . . . . . . . . ∂R l ∂ψ 1 . . . ∂R l ∂ψ l        (3.21)
where system (3.20) is assumed converged to a root {ψ sol } of system (3.19) once {δψ k } → 0. The T mn are the so called tangent matrix terms.

Solution algorithms

Here a brief discussion of three algorithms that can be constructed with a mix of the different physics modules and a choice and/or mix of coupling techniques, is presented in the following:

1. Globally monolithic:

The globally monolithic algorithm refers to an algorithm that solves for a stationary solution by solving for all the modules, defined by their dependent variables ψ = [ϕ, T, B, u] simultaneously.

The nonlinearity of the different modules, requires that the coupled system be solved using the Newton-Raphson technique. The globally monolithic system to solve will be of the form:

   T ϕϕ . . . T ϕu . . . . . . . . . T uϕ . . . T uu    k-1      δϕ δT δB δu      k =        ∂R ϕ ∂ϕ . . . ∂R ϕ ∂u . . . . . . . . . ∂R u ∂ϕ . . . ∂R u ∂u        k-1      δϕ δT δB δu      k = -      R ϕ R T R B R u      k-1 (3.22)
Assuming that the boundary and interfacial conditions are included into system (3.22), solving the linearized system until δψ → 0 is the globally monolithic manner of finding stationary solutions to TIG welding configurations. However, in practice it is difficult and impractical to create the full Jacobian at every iteration k. Moreover, inverting the full Jacobian at every iteration k would be a daunting task for classical solvers, both direct and iterative [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF][START_REF] Viot | Couplage et synchronisation de modeles dans un code scenario d'accidents graves dans les reacteurs nucleaires[END_REF]. Therefore, this algorithm is not used in this thesis.

Partitioned by module:

The modularly partitioned algorithm is one that makes use of the ability of the Newton-Raphson method to linearise nonlinear and implicit equations, while benefiting from the ease of applicability of the Gauss-Seidel technique. This method has a convergence trend that is weaker than the globally monolithic scheme, but is a much more practical algorithm of choice [START_REF] Patricot | Couplages multi-physiques : évaluation des impacts méthodologiques lors de simulations de couplages neutronique/thermique/mécanique[END_REF]. Here, only the diagonal terms of the global Jacobian (3.21) are used, and the linearised, simplified system is partitioned by module into separate linear systems. Finally, the number of linear systems to couple and iterate over using the Gauss-Seidel technique is equal to the number of variables to solve for.

If ψ = [ϕ, T, B, u],
then dim(ψ) = 4 linear systems are solved for k iterations.

                             [T ϕϕ ]{δϕ} k = ∂R ϕ ∂ϕ k-1 {δϕ} k = -{R ϕ (ϕ k-1 , T k-1 , B k-1 , u k-1 )} [T T T ]{δT } k = ∂R T ∂T k-1 {δT } k = -{R T (ϕ k , T k-1 , B k-1 , u k-1 )} [T BB ]{δB} k = ∂R B ∂B k-1 {δB} k = -{R B (ϕ k , T k , B k-1 , u k-1 )} [T uu ]{δu} k = ∂R u ∂u k-1 {δu} k = -{R u (ϕ k , T k , B k , u k-1 )} (3.23)
Assuming that the boundary and interfacial conditions are included into system (3.23), solving the linearised and partitioned system until the set of all variables ψ → 0 is the modularly partitioned manner of finding stationary solutions to TIG welding configurations. In fact, this method is fundamentally a Gauss-Seidel algorithm, and so it is simple to implement, where the coupling occurs by simple data transfer between modules. Furthermore, the fact that the method retains the linearisation capability of the Newton-Raphson method, makes it effective for solving non-linear problems.

The creation and inversion (by direct or iterative methods) of the diagonalised and partitioned tangent matrices is manageable [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF]. Therefore, this algorithm is used to couple the different physics modules used in this thesis. The tangent matrices and the residuals presented in partitioned system (3.23) are detailed in appendix C. Moreover, a physics module can be further partitioned by domain, as is performed for the fluid dynamics module in this thesis.

Partitioned by domain:

The fluid dynamics module treats a two fluid system, where the two fluids are modelled by different hypotheses. Furthermore, the treatment of the fluid interface in this thesis, requires that the fluid dynamics module be solved on a deformable mesh. Additionally, the fluid system can require further partitioning if the linear systems constructed at each iteration k become too large to be practically handled by the linear solvers (by direct or iterative methods). Thus, an algorithm that further partitions a module by domain, is of interest. Partitioning by domain, as it is implemented in this thesis, is treated and discussed in detail in section 3.6. This method also utilises the Gauss-Seidel coupling method, except that the coupling is between domains and not between modules.

The algorithms and coupling techniques that are discussed above are used to solve the global multiphysics model in this thesis. Effectively, both the electric and magnetic modules are solved in a straightforward manner where the linearised system is constructed over the entire TIG domain. As for the conjugate heat transfer and fluid dynamics modules, they are first further manipulated to account for the coupling of the different physical models defined for the arc-plasma and the weld pool. The conjugate heat transfer module accounts for both the electro-thermal phenomena (in the cathode and arc-plasma) and the phase Before proceeding I define the term rigid mesh. It refers to the computational domain before the interface deformation field h z is calculated and applied to the mesh. Due to the partitioned nature of the numerical approach used in this work, the electric, magnetic, temperature and velocity fields are calculated on a rigid mesh, that is later deformed, at every global iteration.

Conjugate heat transfer coupling algorithm

This section discusses the solution methodologies, and the domain/field coupling method used to solve for the heat and phase change phenomena in the TIG system, as implemented in this thesis. The schemes used to solve the energy conservation equations in the Ω cat,arc and Ω pool domains are first presented; followed by a presentation of the coupling scheme used to communicate the temperature approach, used in the arc, to the enthalpic approach, used in the pool, at the fluid interface at Γ AP I .

Interfacial thermal conditions, a mixed heat equation approach

Resolving the heat field of the arc using the standard temperature approach, and with the heat field of the workpiece requiring the enthalpic approach to effectively capture the phase change phenomenon, the arc-pool interface requires careful treatment. The interface between the arc and the workpiece is treated using equation (2.74). In looking to avoid numerical instabilities, the interface is treated by strongly coupling the arc and pool domains [START_REF] Michael | Stability analysis of numerical interface conditions in fluid-structure thermal analysis[END_REF], and the energy conservation equations across the Γ AP I interface are solved monolithically. Before treating the coupling of the interface, the linear systems of the subdomains in question are presented.

Starting with the heat equations in the cathode and the arc-plasma, equation (3.4) is linearised. The linearised system treats the temperature as the primal variable, for iteration k and is expressed in the following:

[T T T ] k-1 {δT } k = -{R T } k-1 (3.24)
where [T T T ] represents the temperature tangent matrix, {δT } and {R T } the temperature increment and residual vectors, defined across the nodes of the meshed cathode and arc subdomains, respectively. As for the heat equation in the workpiece, equation (3.5) is linearised. The linearised system treats the enthalpy as a primal variable, for iteration k and is expressed in the following:

[T hh ] k-1 {δh} k = -{R h } k-1 (3.
T T T 0 0 T hh k-1 δT δh k = - R T R h k-1 (3.26)
The concatenation allows the conjugate heat trasfer module to be solved simultaneously, without the need to transmit temperature or flux fields from one subdomain to another in an iterative manner [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF][START_REF] Michael | Stability analysis of numerical interface conditions in fluid-structure thermal analysis[END_REF]. However, linear system (3.26) is ill-posed in its current form due to the temperature and enthalpy equations being open at their respective Γ plaapi and Γ anoapi subdomains. Thus, with equation (2.74) as a closure condition, and in referring to figure 3.5, the interfacial condition interpreted across the plasma and pool domain is:

T k i = h k j c k p j ∀ i = j (3.27)
and the associated constraint:

λ k i,j = Γ i [(κ∇T ) • n] k dΓ i - Γ j κ∇ h c p • n k dΓ j ∀ i = j (3.28)
where λ k i,j is the Lagrange Multiplier at the overlapping nodes i = j, which manifests as the sum of reaction fluxes at i = j, and n the unit normal at its respective node. Appending system (3.26) with conditions (3.27) and (3.28), the heat transfer scheme to be solved becomes:

[J T,h ] k-1 {δT, h} k =    T T T 0 1 i 0 T hh -1 j 1 i -1 j /c p j 0    k-1    δT δh λ i,j    k = -    R T R h 0    k-1 (3.
δT k i = δh k j c k p j
≈ δT k j which in the limit10 that k → K Global , δT k → ϵ and δh k → ϵ, temperature continutiy at the interface Γ IAP is assumed. Linear system (3.29) is representative of the coupling algorithm used to solve the heat transfer module in the arc-pool system. At every global iteration k, the enthalpy field calculated in the anode is used to calculate the temperature field associated to it. The inverse relation to h = T T ref c p T dT is used to this end, by simply using the inverse of the h-T relations as extracted from the thermophysical data set used in this thesis (see appendix A). This amounts to using the calculated enthalpy fields and the T -h relation, schematised by figure (3.6) to find {T anode }. Figure 3.6 displays a typical T -h relation to be used in calculating the temperature field from the enthalpy, by simple interpolation. This methodology allows for the algorithm to use the fixed mesh at iteration k to identify the pool solidus and liquidus boundaries. The pool domain is then used to solve the momentum and mass conservation equations in the pool subdomain [START_REF] Nedjar | An enthalpy-based finite element method for nonlinear heat problems involving phase change[END_REF][START_REF] Vaughan R Voller | An enthalpy method for convection/diffusion phase change[END_REF]. The identification of the pool subdomain is described in the following subsection.

Identifying the pool boundaries

The mushy zone in the workpiece is of high importance as it is the zone where the latent heat of the melt pool is either captured (during melting) or released (during solidification). With the Stefan number of a typical welding process is of Ste ∼ 1 (see section 1.3), the treatement of the mushy zone can significantly impact the temperature field in the concerned workpiece. A typical finitely sized melt pool is bounded by a solid boundary where the pool velocities approach the velocity of the solid boundary. Without loss of generality, the pool velocity approaches v Γs = 0 m•s -1 as it solidifies in the mushy zone. The no-slip condition applied at the solid boundary implies that the mushy zone functions as a transition zone for the momentum equations [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Vaughan R Voller | An enthalpy method for convection/diffusion phase change[END_REF][START_REF] Kumar | 5.12 -casting simulation methods[END_REF][START_REF] Chien Nguyen | Sensitivity to welding positions and parameters in GTA welding with a 3D multiphysics numerical model[END_REF][START_REF] Traidia | Multiphysics modelling and numerical simulation of GTA weld pools[END_REF]. Therefore, the identification of the pool as a subdomain is set to be a function of the anode temperature field and the velocity field as they are calculated in a simulation. Thus, the pool is numerically defined using the two following criteria:

if Ω pool ⊂ Ω Anode ∀ Ω pool ∪ Ω Anode = Ω Anode then
The temperature based criterion:

∀ x pool ∈ Ω anode and if T pool > T s ⇒ T pool = f (x pool )
The velocity based criterion:

∀ l pool ∈ Ω pool and if ∥v pool ∥ > max(∥v pool ∥) × 10 -4 ⇒ v pool = g(l pool )
Then Ω pool is defined at x pool ∪ l pool = r pool where r pool ∈ Ω pool . Once the pool domain r pool is identified, it can be transferred to the pool momentum and mass transfer module in the global algorithm, and this process is reiterated K Global times. The final subsection in this section presents a proof of concept that shows the numerical stability of the mixed thermal approach.

A note on the stability of the interfacial temperature

The stability of monolithic approaches to the conjugate heat transfer problems for complex systems is well known and the reader is referred to the article by Giles [START_REF] Michael | Stability analysis of numerical interface conditions in fluid-structure thermal analysis[END_REF] for a formal discussion of interfacial stability. In an effort to bypass any possible instabilities related to domain-decomposition or partitioned methods, the arc and anode heat equations are solved monolithically in the globally partitioned multiphysics algorithm of this work. Moreover, the monolithic choice for the conjugate heat transfer module of this work is suitable because the constructed linear systems built at every global iteration k are not as massive as those of a momentum and mass transfer problem, even when solving for 3D problems [START_REF] Sheldon | Methodology for comparing coupling algorithms for fluid-structure interaction problems[END_REF]. Therefore, the monolithic conjugate heat transfer approach in this thesis is treatable by the use of classic linear system solvers (Crout or BICGStab). Linear system (3.29) represents the entirety of the conjugate heat transfer problem solved monolithically in this work; however, the mixed thermal approach used in this thesis is slightly different to the classic approaches tackled in the literature [START_REF] Michael | Stability analysis of numerical interface conditions in fluid-structure thermal analysis[END_REF][START_REF] Viot | Couplage et synchronisation de modeles dans un code scenario d'accidents graves dans les reacteurs nucleaires[END_REF]. Generally, only one primal variable is used to solve a conjugate heat transfer problem as opposed to the two T cat∪pla , h ano used in this work. Furthermore, although a detailed theoretical analysis of the expected stability of the mixed thermal approach is out of the scope of this text, a sample axisymmetric TIG Spot simulation is simulated with this algorithm. The welding parameters of the simulated configuration are, an electric intensity of I = 75 A for a thoriated tungsten cathode, an argon arc and a 316L steel. As seen in figure 3.7.a) the convergence of the temperature profiles (at ≈ 2000 K) at the Γ plaiap and Γ anoiap zones is achieved after only k ≈ 50 global iterations. This indicates that the temperature continuity condition of perfect thermal contact is approached, see equation (2.74). Furthermore, the global algorithm convergence plot is presented in figure 3.7 (b), where all primal fields show a convergent trend. Thus, the reader is provided with a preliminary proof of robustness of the monolithic, mixed thermal approach developed in this thesis.

Momentum and mass transfer coupling algorithm

This section discusses the solution methodologies, and the domain coupling methods implemented in this thesis, used to solve the fluid dynamics module in the arc-pool system. The schemes used to solve the momentum and mass equations in the Ω arc and Ω pool domains are first presented; followed by a presentation of the coupling schemes used to communicate the kinematic and dynamic interactions at the fluid interface at Γ AP I .

The algebraic systems of equations

Focusing on the fluidic effects only in the arc and pool subdomains, the following discussion presents the linearised algebraic systems that numerically treat the momentum and mass transfer equations in their respective subdomains.

The arc scheme

The arc, being modelled as a dilatational fluid, has a linearised system of equations that is different from that of the pool. Linearising equations (3.8) and (3.9), the algebraic system solved at every global iteration k is the following:

[J u,P ] k-1 {δ(u, P )} k = -{R u,P } k-1 (3.30)
where [J u,P ] represents the approximated Jacobian matrix, {δ(u, P )} and {R u,P } the velociy/pressure increment, and the residual vectors respectively, as they are defined across the nodes of the meshed Ω arc subdomain. The chosen scheme is based on direct U -P coupling, solving directly for the velocity and pressure fields simultaneously [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF]. The linear system for the arc fluid equations is comprised of the following scheme, expanding equation (3.30):

T uu T uP T P u 0 k-1 δu δP k = - R u R P k-1 (3.31)
where the details of the different tangent matrices are found in appendix C. The boundary conditions are assumed included in linear system (3.31) and are not detailed for brevity.

The pool scheme

The weld pool being modelled as an incompressible fluid has a linearised system different than that of the arc. Linearising equations (3.10) and (3.11), the algebraic system solved at every global iteration k is the following:

[J v,p ] k-1 {δ(v, p)} k = -{R v,p } k-1 (3.32)
where [J v,p ] represents the approximated Jacobian matrix, {δ(v, p)} and {R v,p } the velociy/pressure increment, and the residual vectors respectively, as they are defined across the nodes of the meshed Ω pool subdomain. The chosen scheme is also based on direct U -P coupling [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF]. The linear system for the pool fluid dynamics equations is comprised of the following scheme, expanding linear system (3.32):

T vv T vp T pv 0 k-1 δv δp k = - R v R p k-1 (3.33)
where the details of the different tangent matrices are found in appendix C. The boundary conditions are assumed included in linear system (3.33) and are not detailed for brevity.

Coupling schemes at the fluid interface

The fluid interface at the Γ AP I subdomain is a zone where the kinematic and dynamic conditions are imposed. These conditions allow for the transfer of physical quantities between the arc and the pool and are necessary for the mathematical coupling of the algebraic systems (3.31) and (3.33). The two coupling methods studied in this work are the Dirichlet-Neumann method 11 and a quasi-monolithic method [START_REF] Sheldon | Methodology for comparing coupling algorithms for fluid-structure interaction problems[END_REF].

The coupling methods are applicable to the two interface hypotheses studied in this thesis, the fixed and deformable interface hypotheses. In the following discussion, the details of the coupling techniques for a rigid mesh calculation at every global iteration k are presented first, and then followed by a discussion of the treatment of the fluid interface deformation.

The Dirichlet-Neumman algorithm

The Dirichlet-Neumann coupling algorithm is a classic method, adapted from the literature on numerical fluid structure interaction and conjugate heat transfer [START_REF] Sheldon | Methodology for comparing coupling algorithms for fluid-structure interaction problems[END_REF][START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF][START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF][START_REF] Michael | Stability analysis of numerical interface conditions in fluid-structure thermal analysis[END_REF]. The method consists of decomposing the fluid domain into two partitioned subdomains Ω Arc ∪ Ω pool = Ω f luid , separated by the interface Γ AP I ;

and where each partition Ω Arc and Ω pool is associated to a separate calculation step. Thus, the coupling scheme consists of a total of two steps, where information communicated between the arc and pool cal- 1. The arc fluid system is appended, at its interfacial nodes, by mapping the velocity values from iteration k -1 from Γ anoapi to Γ plaapi at iteration k; and applying the mapped velocities like Dirichlet boundary conditions using the Lagrange Multiplier technique [START_REF] Pellet | Dualisation des conditions aux limites[END_REF]. Figure 3.9 schematises the applied velocity mapping, with the red arrows representing the mapped Dirichlet conditions. T uu T uP I 3,i T P u 0 0

∀ i, j ∈ [1, N ] ∧ ∀ i = j ⇒ ∃ r j ∈ Γ anoapi ⇒ u k i = v(r j ) k-1 (3.34) α k i = Γ i [-P + µ(∇u + ∇ T u - 2 3 ∇ • uI) • n] k dΓ i (3.
I 3,i 0 0    k-1    δu δP α i    k = -    R u R P u i    k-1 (3.36)
2. Once linear system (3.36) is solved, the calculated forces from the arc are translated and transmitted to the pool fluid dynamics scheme, using the Lagrange Multipliers α k i that were solved for. Figure 3.9 represents a visualisation of the applied reaction force mapping, schematised by the blue arrows, representing the transmitted forces. The translation of the fields is applied numerically in the computer algorithm, where the fields are translated from the Γ plaapi nodes to the Γ anoapi nodes. Where for the rigid pool system 12 , the mapped and transmitted arc forces and the impermeability condition, from equation (2.43), at the interface are:

∀ i, j ∈ [1, N ] ∧ ∀ i = j ⇒ ∃ r i ∈ Γ plaapi ⇒ F k τ,j = -α(r i ) k • τ k j (3.37) v k j • n k j = 0 (3.38)
β k n j = n k j • Γ j (n • [-p + µ(∇v + ∇ T v) • n]) k dΓ j (3.39)
where n k i and τ k i are the unit normal and tangent vectors at node i at iteration k. Applying the forces from equation (3.37) to linear system (3.33) similar to Neumann conditions; and appending linear system (3.33) with equations (3.38) and (3.39); the second linear system to solve in this coupling algorithm becomes: m is limited to a maximum of 3 iterations, to help limit redundant calculations, and move on to the next physics modules in the global algorithm. Thus, solving for stationary solutions, the limit is set so that as k → K Global , δu k → ϵ and v k → ϵ, satisfies and ensures a stationary solution, and the kinematic and dynamic conditions at the interface Γ IAP .

   T vv T vp n j T pv 0 0 n j 0 0    k-1    δv δp β n j    k = -    R k-1 v + F k τ,j R k-1 p 0 j    (

A note on the stability of the interfacial velocities

The presented Dirichlet-Neumann scheme has been tested for stability for a wide range of welding configurations 13 , and the velocities tend to converge to the no-slip conditions after a couple of iterations. Figure 3.11 presents the maximum velocity at Γ plaapi and Γ anoapi as they converge rapidly with respect to the algorithm global iteration number k. The displayed results are extracted from the same sample simulation discussed in section 3.5.3. The convergence of the velocity profiles (at ≈ 0.175 m•s -1 ) at both the Γ plaapi 13 Details not discussed here. and Γ anoapi subdomains is acheived quite rapidly, with their error dropping below 1 % in the first couple of iterations. This provides the reader with a preliminary proof of the effectiveness of the Dirichlet-Neumann coupling scheme in solving the momentum and mass transfer problem presented in this thesis.

The quasi-monolithic algorithm

The quasi-monolithic coupling scheme is a more natural14 scheme to solving the two-fluid problem. This approach strongly imposes the interfacial conditions at the Γ AP I interface. The scheme enforces the noslip, see equation (2.44), and impermeability conditions at every iteration k, unlike the Dirichlet-Neumann scheme that imposes them weakly, and that can require subcycling to satisfy the physical interfacial conditions. To apply the no-slip, and the arc and pool impermeability conditions numerically, the Lagrange Multiplier technique is used [START_REF] Pellet | Dualisation des conditions aux limites[END_REF], in:

No-slip:

u k i -v k j = 0 (3.41) λ k i,j = Γ i [-P + µ(∇u + ∇ T u - 2 3 ∇ • uI) • n] k dΓ i - Γ j [-p + µ(∇v + ∇ T v) • n] k dΓ j (3.42)
Impermeability:

v k j • n k j = 0 (3.43) ζ k n j = n k j • Γ j (n • [-p + µ(∇v + ∇ T v) • n]) k dΓ j (3.44)
Assembling the monolithic system, first, linear systems (3.31) and (3.33) are concatenated.

    
T uu T uP 0 0

T P u 0 0 0 0 0 T vv T pv 0 0 T vp T pp      k-1      δu δP δv δp      k = -      R u R P R v R p      k-1 (3.45)
The concatenation allows the momentum and mass conservation laws of the entire rigid fluid domain to be solved simultaneously, without the need to transmit velocity and force fields from one subdomain to another [START_REF] Sheldon | Methodology for comparing coupling algorithms for fluid-structure interaction problems[END_REF][START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF][START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF][START_REF] Viot | Couplage et synchronisation de modeles dans un code scenario d'accidents graves dans les reacteurs nucleaires[END_REF]. Closing the system of equations in system (3.45), equations (3.41), (3.42), (3.43), (3.44) are appended to it; where the monolithic scheme for the arc-pool fluid system for a rigid interface becomes, at iteration k:

         T uu T uP 0 0 I 3,i 0
T P u 0 0 0 0 0 0 0 T vv T pv -I 3,j n j 0 0 T vp T pp 0 0 I 3,i 0 -I 3,j 0 0 0 0 0 n j 0 0 0          k-1          δu δP δv δp λ i,j ζ n j          k = -          R u R P R v R p 0 i,j 0 j          k-1 (3.46)
Linear system (3.46) represents the monolithic scheme as it is applied at at every iteration for the rigid interface, fluid dynamics module. Thus, solving for stationary solutions, the limit is set so that as k → K Global , δu k → ϵ and v k → ϵ, satisfies and ensures a stationary solution.

A note on stability:

As per the literature, monolithic methods have no instabilities associated with interface treament. This makes the stability of the spatial discretisation of the system both necessary and sufficient in ensuring the numerical stability of the algorithm. This is unlike Dirichlet-Neumann coupling, where stability depends on spatial discretization, the coupling direction and other criteria [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF][START_REF] Michael | Stability analysis of numerical interface conditions in fluid-structure thermal analysis[END_REF]. The reader might thus assume that monolithic methods are to be preferred, however, monolithic methods can easily grow in degrees of freedom and can rapidly exhaust linear solvers [START_REF] Viot | Couplage et synchronisation de modeles dans un code scenario d'accidents graves dans les reacteurs nucleaires[END_REF][START_REF] Sheldon | Methodology for comparing coupling algorithms for fluid-structure interaction problems[END_REF][START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF]. Partitioned methods like the implemented Dirichlet-Neumann method are thus useful, but should be tested for stability thoroughly before use.

A note on the rigid interface and quasi-monolithicity:

The discussion up until now has only presented the treatment of the fluid dynamics problem as it occurs for a rigid mesh. That is, both the Dirichlet-Neumann and the quasi-monolithic algorithms presented in this section couple the arc and pool domains without considering the deformation at the interface. This, motivates the presentation of the following brief discussions:

• Simulating welding using the fixed interface hypothesis, the unit normals n used in equations (3.40) and (3.46) calculated at every iteration k15 need only be assigned the initial rigid interface orientation, n initial = (n x,0 , n y,0 , n z,0 ). No treatment of the rigid interface dynamics, nor deformation vector is thus needed.

• The coupling algorithms that have been presented until now implement the coupling at the rigid Γ AP I interface. This is because the treatment of the deformable interface problem in this thesis is performed in a partitioned manner. This brings into context the notion of quasi-monolithicity when working with deformable interfaces [START_REF] Sheldon | Methodology for comparing coupling algorithms for fluid-structure interaction problems[END_REF]. Effectively, the calculation on the rigid fluid mesh is performed monolithically, where the interfacial forces are then transferred to the interface deformation calculation.

Fluid interface deformation

The way the fluid interface at Γ AP I is treated in this thesis is based on the trial method 16 , which is essentially a partitioned method [START_REF] Cuvelier | Some numerical methods for the computation of capillary free boundaries governed by the navier-stokes equations[END_REF]. Being a partitioned method, this approach is compatible with the coupling schemes implemented in this thesis. The trial method consists of first solving for the Navier-Stokes equations on a rigid mesh, and then solving for the interface deformation field h z , by looking for its equilibrium state. Once h z is calculated it is used to deform the arc and pool meshes in a regularised manner. It is noteworthy to state that the interface deformation is calculated at Γ anoapi nodes and then mapped to the Γ plaapi nodes (see figures 3.8, 3.9). The following cases discuss the arc and pool reaction forces, as per the coupling scheme is used:

• Dirichlet-Neumann:

If linear systems (3.36) and (3.40) are solved for, the translation of the arc forces at Γ plaapi to Γ anoapi is done in the same way represented by figure 3.9. The arc and pool forces at node j becomes:

F k arc,n,j = -α(r i ) k • n k j F k pool,n,j = -β k n j
• Quasi-monolithic:

If linear system (3.46) is solved for, the translation of the arc forces at Γ plaapi to Γ anoapi is done in the same way represented by figure 3.9, the arc and pool forces at node j becomes:

F k arc,n,j = -(λ(r j ) k • n k j ) + ζ k n j F k pool,n,j = -ζ k n j
where once F arc,n,j and F pool,n,j are calculated, they are added to the sum of normal forces at Γ AP I , in:

R dyn,j = (F arc,n,j -F pool,n,j -F σ,j -F g,j ) • ẑ (3.47)
where R dyn,j is the dynamic residual vector defined at j, projected along the vertical axis 17 . This projects the deformation vector x to the vertical axis, making the deformation field to calculate for the h z field [START_REF] Chien Nguyen | Sensitivity to welding positions and parameters in GTA welding with a 3D multiphysics numerical model[END_REF][START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF]. Furthermore, to ensure the closure of the dynamic system of equations, the volume conservation equation at the interface, imposed similarly to a Dirichlet boundary condition, is numerically formulated as:

R kin = L l=1 V l -V 0 (3.48) V l = 1 dim(R n ) Γ l x • n N hz dΓ l (3.49)
where R kin is the kinematic residual vector, V l the volume of element l for L elements at the interface, and

where

V 0 = L l=1 V l k=0
, being the initial pool mesh configuration volume. The solution strategy used to solve the interface deformation problem is the following:

1. The interface deformation problem is linearised, as before, using the quasi-Newton method, to give the following system of equations:

T hzhz M T hzhz M hzhz 0 k-1 δh z ψ k = - R dyn R kin k-1 (3.50)
where the details of the tangent matrices are presented in C.

2. After solving for the deformation field h k z,j , defined at the Γ anoapi subdomain, the field is transmitted to the geometrically intersecting nodes at i of the Γ plaapi subdomain. Referring to figure 3.9, the deformation field is copied to h k z,i , in:

∀ i, j ∈ [1, N ] ∧ ∀ i = j ⇒ ∃ r j ∈ Γ anoapi ⇒ h k z,i = h z (r j ) k (3.51)
3. With the calculated deformation vector loaded at the interface nodes, the arc and pool internal elements, at Ω arc and Ω pool , need displacement, following the deformation of the interface Γ AP I . To this end, a mesh regularization method is used, as adapted from Huang [START_REF] Huang | Mathematical principles of anisotropic mesh adaptation[END_REF]. The basic principles of this method are presented in the following:

For every internal element18 Ω m , an element energy is associated to it, described by:

E m = E m1 + E m2 (3.52)
where:

E m1 = θ H 2 Ωm tr ([G T ][M ][G]) -1 nγ H 2 det[M ] dΩ m E m2 = (1 -θ H )n nγ H 2 Ωm det[G T ][M ][G] 1-γ H 2 dΩ m
where θ H is a weighting factor chosen to be 0.5, and [G] is the Jacobian of the geometric transformation between element Ω k-1 m and a reference element Ω k m at iteration k, and [M ] is the target metric tensor, and n = dim(R n ), and γ H = 2 representing the method norm. The mesh displacement method used in this thesis does not look to adapt the mesh, rather, only to regularize it w.r.t the interface deformation field {h z }. Thus, the target metric tensor is:

[M ] = ([G 0 ][G T 0 ]) -1 where [G 0 ] = [G] k=0
which constrains all element deformations at iteration k to be similar to the initial rigid element at the start of the calculation at k = 0. The first element energy contribution E m1 19 of equation (3.52) controls the regularity of the element, while the second contribution E m2 controls the volume of the element. The internal element displacement fields are calculated at every iteration k, by: a) Calculating the node displacement field vector {δX pool } k of the pool domain by setting, and constraining the system with h z,j at the Γ anoapi subdomain, for E m ∈ Ω pool :

   ∂ 2 E m ∂X 2 pool 1 j 1 j 0    k-1 δX pool ϑ j k = -   ∂E m ∂X pool δh z,j   k-1 (3.53)
once {δX pool } k is found, the pool nodes are displaced by a simple mesh position update

{X pool } k = {X pool } k-1 + {δX pool } k b)
Calculating the node displacement field vector {δX arc } k of the pool domain by setting, and constraining the system with h z,i at the Γ plaapi subdomain, for E m ∈ Ω arc :

   ∂ 2 E m ∂X 2 arc 1 i 1 i 0    k-1 δX arc ϑ i k = -   ∂E m ∂X arc δh z,i   k-1 (3.54)
once {δX arc } k is found, the arc nodes are displaced by a simple mesh position update

{X arc } k = {X arc } k-1 + {δX arc } k
The presented mesh deformation and regularisation method discussed is generally adapted to the partitioned schemes used in this work. Moreover, the method strictly conserves volume throughout the interface deformation calculation, because of the mass conservation constraint (3.48) imposed onto linear system (3.36). A representation of interface deformation, node displacement and regularization algorithm is presented in figure 3.12. Figure 3.12 (a) presents the rigid mesh where a 2D axisymmetric spot calculation is performed, whereas (b) presents the deformed mesh as it is influenced by the fluid interface dynamics. The displayed meshes are extracted from the same sample simulation discussed in section 3.5.3 and 3.6.2.1, but with an artificially magnified arc pressure (F arc,n × 15) for visualisation purposes.

Convergence criteria

Seeing as all the algorithms used to approximate solutions to the non-linear equations on the finitely sized meshes of this work are limited to their discrete nature, finitely sized convergence criteria are necessary. In multiphysics problems, an approximate solution is considered converged once the residuals or the increments of the key primal variables become sufficiently small. The definition of "sufficiently small" depends on the application and on the problem requirements. In general the sufficiently small criterion is based on the idea that a solution is considered converged once the calculated approximate solutions no longer vary per iteration in a simulation. The convergence criteria used in this thesis are based on the increment size of the variable in question. Furthermore, a different set of convergence criteria are defined for the 2D and 3D simulations of this work. This criteria for 3D simulations are less stringent so as to reduce calculation time. Therefore, to be able to monitor the convergence trends of a simulation, a set of key variables are used to calculate the increment at every global iteration. The increments of the voltage, temperature, arc and pool velocities and the interface deformation variables are calculated using the following:

δϕ k = ϕ k -ϕ k-1 |ϕ k max | ; δT k = T k -T k-1 T k max ; δu k = u k -u k-1 |u k max | ; δv k = v k -v k-1 |v k max | ; δh k z = h k z -h k-1 z |h k z max | . (3.55)
The increments seen in the set of equations (3.55) are calculated at every global iteration k. The increments are compared at every k to the defined convergence criteria. The defined convergence criteria for 2D and 3D simulations are presented in the table below: These criteria are used throughout the entirety of this thesis unless otherwise stated explicitly. As a final note, a convergence criterion for the magnetic field module is not used because of the linearity of the magnetic system of equations. This means that during a simulation, the magnetic field increment δB k generally follows the convergence trend of the voltage increment δϕ k .

δϕ min δT min δu min δv min δh z min 2D 1 × 10 -4 1 × 10 -4 1 × 10 -4 1 × 10 -3 1 × 10 -4 3D 1 × 10 -4 1 × 10 -3 1 × 10 -3 1 × 10 -3 1 × 10 -3

Introduction

In this chapter, the study and analysis of the sensitivity of a TIG arc to the emitting cathode geometry is studied and discussed. A 2D axisymmetric TIG Spot configuration is set up and the cathode tip parametrised, and the thermal and dynamic responses of the arc analysed. Although cathode tip geometry is known to influence arc-plasmas, a thorough quantification of the Lorentz phenomenon, which is the dominant source of momentum in the fluid flow of TIG arcs has not been strongly investigated. Therefore, a numerical parametric study is performed at a constant inlet electric current and arc height, for different cathode sizes and shapes (pointed, chamfered and rounded tips), parameterised by the truncation angle and tip radius. A characterisation, quantification and an analysis of the influence of tip geometry on TIG arcs and the workpiece is presented. However, to reduce computational costs and to isolate the influence of the cathode onto the arc, the workpiece is modelled as a solid copper anode.

Geometric configuration

This section presents the problem geometry for a TIG Spot configuration and its parametrisation, assuming a rotational symmetry, as sketched in figure 4.1. The model accounts for a cathode, arc plasma and anode, in which both electrodes are assumed to remain in their solid state to highlight the behaviour of the arc plasma as a function of only the cathode geometry. So, the arc height is set to a constant value of 10 mm throughout the study. The parametric study is concerned with the effect of the variation of the cathode geometry at its tip. The considered shapes are a pointed cathode (PNT), a chamfered one (CHF) and a rounded (RND) one. The cathode is modelled as solid thoriated tungsten (2% Th), the arc as an argon plasma and finally the anode as a solid copper piece. The geometry of the cathode is parameterised

h cat1 + h cat2 = 10 mm r int = [0.15 mm, 0.3 mm, 0.6 mm] α = [15 o , 20 o , 25 o , 30 o ] β α = 90 -α h cat1 = 2 -r int tan α (4.1)
for all cases

r int,rnd = r int tan(β α /2) (4.2) rounded case r int,chf = r int (4.3) chamfered case

Boundary conditions

The imposed boundary conditions for the partial differential equations used in this study are presented in where, the applied inlet electric current j imp at segment [EF] is imposed in the following manner:

Material properties

The thermophysical and electric properties of the solid tungsten cathode, the argon plasma and the solid copper anode are assumed to depend only on temperature. The thermophysical properties (ρ, µ, λ and c p ) and the electrical one (σ * ) are drawn from Brochard's PhD [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF], and reported in appendix A for self-consistency of the present paper.

Discretisation and meshing

This subsection discusses the chosen discretisations used for the different primal variables, the different meshes used and the spatial convergence of some key variables that the arc model solves for. The discussed mesh dependent variable of interest is the arc velocity. Once an appropriate mesh is identified, it is used for similar configurations throughout the rest of this section.

Electric ϕ = 0 [HI] -σ∇ϕ • n = 0 [EA] ∪ [AB] ∪ [BH] ∪ [IC] ∪ [CD] ∪ [DF] -σ∇ϕ • n = j imp [EF] Heat transfer T = 1000 K [EF] ∪ [FD] ∪ [CI] ∪ [IH] -λ∇T • n = 0 [EA] ∪ [AB] ∪ [BH] ∪ [DC] Magnetic B = 0 [EA] ∪ [AB] ∪ [BH] Momentum u • r = 0 [AB] ∪ [BC] u • ẑ = 0 [BC] ∪ [CD] (σ • n) • r = 0 [FD] ∪ [CD] (σ • n) • ẑ = 0 [FD] ∪ [AB]

Discretisation

The primal variables (ϕ, T, B, u) are discretised using (Q2, P2) elements, while the temperature depenendent physical parameters (σ, λ etc...) and the pressure variable P are discretised using (Q1, P1) elements.

The choices of element discretisation per physical variable and coefficient are presented and summarized in table 4.

The different meshes

Spatial convergence

A formal discussion on spatial error is out of the scope of the present study, owing to the use of unstructured meshes [START_REF] Christopher | Grid convergence error analysis for mixed-order numerical schemes[END_REF]. However, I have estimated the spatial discretisation error thanks to 4.4:

relative difference = 1 - f M 2 f M 3 • 100 (4.

Discussion and results

This section presents the results of the parametric study and their analysis for an arc height of 10 mm and a input electric current of 200 A. I discuss from a physical point of view, the obtained results (influence of electric intensity on transported quantites, arc pressure and shear, heat flux).

Physical discussion

In this parametric study aims the magnetic pinching effects induced by the studied geometries of the cathode tip are highlighted. The cathode, arc plasma and the anode are all taken into account in the present model. However, both electrodes are assumed to remain in their solid state. The simplified electric model at the electrode-plasma interfaces allows for a reduction of the computational cost, meanwhile achieving a satisfactory physical representativity. This enables us to account for the geometric effects of the cathode onto the induced electric current density. The arc height of the configuration has been chosen such that the electrical effects at the API become invariant with respect to the studied cathode geometries.

Furthermore, to achieve a good trade-off between an acceptable physical representativity and an affordable computational cost, we have used a purely conductive model in the anode. This aims to model the thermal and electrical sinks in a very cheap way. However this modelling choice infers that any reference to the thermal field in the anode is unrealistic.

Influence of cathode shape on arc behaviour

The Lorentz force f Lor is the leading momentum source term in TIG arcs [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Hsu | Study of the free-burning high-intensity argon arc[END_REF][START_REF] Tanaka | A unified numerical modeling of stationary tungsten-inert-gas welding process[END_REF] which the intense Lorentz force acts on the arc is also much higher for the pointed case. This change in the Lorentz force accordingly influences the velocity field magnitude in the arc-plasma as is seen in the lower row of figure 4.5. In fact, the magnitude of the velocity field rises 2.8 times from the rounded to the pointed cases presented. This results from the intense magnetic pinching effect of the arc at the sharply pointed tip, which leads to a stronger acceleration of the arc in the direct vicinity of the cathode tip. Therefore, the sharper the cathode tip, the higher the constriction of the arc due to increasing Lorentz forces. These effects translate into higher transported heat flux values and imparted stresses at the anode. We note that the considered pointed cathode configurations are ideal geometries and are physically unrealistic. This is due to geometrical singularities at the tip [START_REF] David | Classical electrodynamics[END_REF], and this results in temperatures beyond the melting point of tungsten (T melting = 3680 K) in its vicinity. This is represented in figure 4.6 where the temperature profiles at the cathode symmetry axis for the sharpest and bluntest cathodes are plotted. However, although the pointed cathodes are physically unrealistic, they correspond to an extreme geometric limit of cathode shapes. Being the upper limit of cathode sharpness, the pointed cases set the bound on the behaviour of the implemented physical model; while the bluntest chamfered and rounded cases set the lower bound. Thus, studying the physical model at the geometric extremes, the behaviour of the model at any intermediate cathode geometry can be understood.

Influence of the cathode shape at the plasma-anode interface

To better understand the underlying mechanisms at play in this parametric study, I discuss how the heat flux and mechanical stresses act at the API. Three heat transfer modes make up the total heat transferred to the anode at the API, and they are made up of a diffusion, radiation and an electronic absorption term.

Their computed profiles along the radial direction at the API are presented in figure 4.7. Integrating the corresponding fluxes over the API, the radiated heat power is about 65 W while both the heat powers by electronic absorption and the conducted heat from the plasma are of the order of 10 3 W. These trends are found for all studied cases and not just for the presented case. This implies that only the electronic absorption and conducted heat terms play a dominant role in contributing to the heat flow at the API.

Therefore, the total heat flow absorbed by the anode is dominated by both the electronic absorption and conducted heat terms. The electronic absorption term, being directly proportional to the electric current density at the API, varies proportionally to the j AP I profile. Furthermore, because the conducted heat at the API is influenced by the advection of heat in the arc, it increases as the advection in the arc rises. This The color key groups cathodes of same angle, while the marker key groups cathode shape and tip size.

The rise in conducted heat as a function of arc Peclet number is essentially a rise in the transported heat of the arc as a function of cathode tip sharpness.

In the following subsection, I have quantified the influence of the cathode geometry on to the physical quantities of interest at the API (arc pressure, arc shear stress and the total heat flux to the anode).

Transported quantities at the plasma-anode interface

The impact of the variability in the Lorentz forcing of the flow is directly seen in the arc velocity profiles, where the velocity fields get stronger as the cathode tips get sharper. The direct influence of cathode sharpness on arc velocity is presented in figure 4.9, where the maximum arc velocities are plotted as functions of truncation angle and tip size. This significant rise in maximum arc velocity with tip sharpness significantly influences the pressure, shear and heat flux profiles of the arc at the API. To quantify the variability in the stresses and heat flux profiles at the API, two analytical functions are used to simplify the discussion. The functions are fit to the profiles of interest at the API similarly to Alvarez et al. [START_REF] Delgado-Álvarez | Dimensionless representation of the column characteristics and weld pool interactions for a DC argon arc[END_REF]. However, I use simpler (less parameters) analytical fits, motivated by the work in [START_REF] Kim | Modeling of temperature field and solidified surface profile during gas-metal arc fillet welding[END_REF][START_REF] Sung | Mathematical modeling of the dynamic behavior of gas tungsten arc weld pools[END_REF][START_REF] Chien Nguyen | Sensitivity to welding positions and parameters in GTA welding with a 3D multiphysics numerical model[END_REF].

Truncation Angle

The heat flux q ano and the pressure P at the API are fitted with the following Gaussian distribution:

g(r) = g max exp -r 2 2σ g 2 (4.5)
The shear stress τ at the API is fitted with a scaled Weibull distribution as follows:

w(r) = w max r λ w k k -1 1 k k-1 exp k -1 k - r λ w k (4.6)
The functions are fit by centering them at the maximum values of the profiles q max , P max , τ max and by constraining the integrals of the fitting functions to the integral heat and force values. The widths of the numerical profiles are quantified using the variance σ g of equation 4.5 and the scale λ w of equation 4.6.

The shape parameter k of equation 4.6 dictates the shape of the Weibull distribution.

The variability in the magnetic pinching effect as a function of tip sharpness translates to significant variability in the transported arc quantities at the API. Between the bluntest and the sharpest cathodes, the maximum pressure P max and shear stress τ max at the API increases ≈ 3.5 and ≈ 4-fold, respectively, while the maximum heat flux q max increases ≈ 1.7-fold. Moreover, as the maxima rise with respect to tip sharpness, the widths of the pressure, shear and heat flux profiles consequently drop. The widths σ P , λ τ of the pressure and shear profiles at the API drop ≈ 1.6-fold, while the widths σ q of the heat flux profiles drop ≈ 1.3-fold between the bluntest and sharpest cathode tip. As the stress profiles intensify with cathode tip sharpness, the total imparted integrated pressure and shear forces increase ≈ 1.5-fold and ≈ 1. The sharper the geometry of the cathode tip, the more energy and momentum is transported by the arc from the cathode to the anode. These trends are captured by our model because the cathode is included in the calculation domain. Thus the geometry of a cathode requires attention when setting up a TIG welding configuration for simulation.

Analysis of the obtained results

This study assumes a solid anode so as to isolate the effects of cathode tip geometry onto TIG arcs. The use of a solid anode maintains a constant arc height unlike in actual welding conditions which leads to a deformable melt pool in the anode domain [START_REF] Chien Nguyen | Sensitivity to welding positions and parameters in GTA welding with a 3D multiphysics numerical model[END_REF]. By solely increasing the sharpness of a cathode tip, the momentum and heat transported by a TIG arc increases. Referring to the figures 4.10 and 4.11 I note that the maximum stresses and the heat fluxes at the API increase roughly 4 and 1.7 times, respectively, from the bluntest to the sharpest cathode tips. This variability is significant and important to capture in TIG welding simulations, as is already reported in the literature on experimental observations [START_REF] Key | Anode/cathode geometry and shielding gas interrelationships in gtaw[END_REF][START_REF] Savage | The effect of electrode geometry in gas tungsten-arc welding[END_REF][START_REF] Mills | Factors affecting variable weld penetration[END_REF][START_REF] Alber | Gas-tungsten-arc cathode and related phenomena[END_REF][START_REF] Tsai | Heat distribution and weld bead geometry in arc welding[END_REF].

Therefore, when setting up numerical simulations, one should pay attention when accounting for cathode shape in the model. Indeed, it appears that in some early numerical studies the sensitivity of the momentum and mass equations to the magnetic pinching effect is overlooked [START_REF] Hsu | Study of the free-burning high-intensity argon arc[END_REF][START_REF] Goodarzi | The effect of the cathode tip angle on the GTAW arc and weld pool: I. mathematical model of the arc[END_REF][START_REF] Bini | Numerical and experimental study of transferred arcs in argon[END_REF]. With the Lorentz force being the main source of momentum in the arc-plasma, any variation in this force translates into a variation in the velocity field, as is presented in figure 4.9. Furthermore, the variability of the Lorentz force is directly due to its sensitivity to the electric current density in the vicinity of the cathode tip. From figure 4.4, I note the sensitivity of the electric current density to the geometry at the cathode tip. This sensitivity is captured by the explicit inclusion of the cathode into the calculation domain. The sensitivity is due to the ∇ϕ| Γ CP I term that increases in intensity as the tip geometry gets sharper [START_REF] Fix | An analysis of the finite element method[END_REF][START_REF] David | Classical electrodynamics[END_REF].

For the considered arc-plasma of a 10 mm height and 200 A inlet current, the calculated variability in the heat flux transported to the API is dictated by the variability in the conducted heat term q pla-ano . This is because the variability in the electronic absorption term s Abs at the API is negligible due to the arc becoming relatively insensitive to the voltage field near the anode. This insensitivity implies that the electric current density once it reaches the API becomes roughly constant, as shown in the inequalities (4.7) and (4.8). The calculated conducted heat flux q pla-ano varies with the sharpness of the cathode tip due to its dependence on the advected heat of the arc. Referring to figure 4.8, it turns out that the transferred heat flux to the anode depends on the Peclet in a quasi-linear way. Therefore, the magnitude of the arc velocity, which depends on the Lorentz force and consequently on the cathode geometry, strongly influences the heat flux at the API.

Although variability is calculated in the transported quantities across most studied cases, negligible variation is seen between simulations of the CHF and RND type; as is observed in figures 4.8, 4.9, 4.10, Furthermore, similar results are calculated for multiple cases at different truncation angles and radii. This is observed in figures 4.8, 4.9, 4.10, 4.11. This weak variability between the cases is attributed to a quasi-equivalence in the total power of the Lorentz force to drive flow in the arc-plasma. The total power is defined by:

P = Ω pla u • f Lor dΩ (4.9)
where P [W] represents the total power in driving the fluid flow. The power of the Lorentz force for all simulated arcs in this study is presented in figure 4.13. Referring to the figure, the cases with quasi-equivalent power values are encased by a specific black marker (circular, triangular and rectangular markers). Provided that the integral power of the Lorentz force is similar, the resulting transported quantities at the API will also be comparable, irrespective of the cathode shape.

Voltage drop and arc efficiency

The voltage drop across the TIG circuit is calculated by imposing a constant current at the input boundary of the cathode.

η = Q ano P circuit + s Abs -s Emi + s N eu dΓ × 100% (4.10)

A discussion with respect to similar works

In this section the results of the present study are compared to those found in the literature on both numerical and experimental works.

Comparison to simulations

The first comparison is performed with a numerical study from Goodarzi et al. [START_REF] Goodarzi | The effect of the cathode tip angle on the GTAW arc and weld pool: I. mathematical model of the arc[END_REF], in which a wide range of truncation angles 1 , ranging from 4.59 o ≤ α ≤ 75 o , were considered. Instead of incorporating the cathode into their calculation domain, they imposed an electric current density at the cathode tip as a boundary condition. The electric current density is imposed using spot areas experimentally estimated by Haidar et al. [32]. translates into negligible differences in the transported quantities of the arc. For ease of comparison and due to the fact that the spot size is the critical boundary condition in Goodarzi's work [START_REF] Goodarzi | The effect of the cathode tip angle on the GTAW arc and weld pool: I. mathematical model of the arc[END_REF], the cathode spot area sizes of this study are estimated. I define the spot size as the zone where 150 A of current passes through the cathode. The estimated spot areas are presented in figure 4.15 figure (b), along with those from Goodarzi [START_REF] Goodarzi | The effect of the cathode tip angle on the GTAW arc and weld pool: I. mathematical model of the arc[END_REF] and Haidar [32]. The variability of the spot sizes w.r.t truncation angle, tip radius and shape shows a ≈ 3 fold maximum variation in area size. Moreover, the estimated spot areas do not exceed 2.4-times the values used by Goodarzi et al.. Thus, the variability calculated in the estimated spot areas is not sufficient to explain the ≈ 1000-fold (see table 4.4) ratio that I obtained in the electric current densities between the bluntest and sharpest cathode tips. Therefore, I deduce that the spot area boundary condition method used by Goodarzi et al. is not adapted for the capturing of the variation of electric current density at the tip. Indeed, the spot areas are incapable of capturing the local behavior at the cathode tip, because the spot is interpretable as an integral quantity; masking any local j gradients that can influence the Lorentz force of the arc flow. Thus, the spot area boundary condition method is not adapted to capturing the influence of tip sharpness onto TIG arcs. Thus, the necessity to account for the cathode geometry in the computational domain when modelling TIG arcs is highlighted.

Variable [START_REF] Goodarzi | The effect of the cathode tip angle on the GTAW arc and weld pool: I. mathematical model of the arc[END_REF] and from the present study.

j CP I,max [A•mm -2 ] q max [W•mm -2 ] u max [m•s -1 ]

Comparison to experiments

In this subsection, I compare the calculated results to several experimental studies that deal with the influence of cathode tip sharpness on to TIG arcs. Mills et al. [START_REF] Mills | Factors affecting variable weld penetration[END_REF] note the 4-fold increase in maximum arc pressure as the cathode tip truncation radius decreases to r int = 0 for a constant tip truncation angle of α = 15 o . Furthermore, the maximum arc shear stress rises similarly with the arc pressure for typical TIG arcs. Tsai [START_REF] Tsai | Heat distribution and weld bead geometry in arc welding[END_REF] observes that the maximum heat flux transported to the anode increases 1. [START_REF] Mills | Factors affecting variable weld penetration[END_REF] ≈ 390

≈ 1470 ≈ 4 -- -- ≈ 4 -- -- -- Tsai [86] -- -- -- -- -- -- ≈ 27 ≈ 32 ≈ 1.2 Petrie [87] ≈ 400 ≈ 1400 ≈ 3.5 -- -- -- -- -- 1.3 Sadek [71] ≈ 500 ≈ 2200 ≈ 4.5 -- -- -- -- -- --

Summary

The numerical parametric study showed the influence of the cathode tip geometry onto relevant welding quantities, such as the heat transferred to the anode, the arc pressure and shear stress at the workpiece.

The TIG Spot model implemented in the present work accounts for both the cathode, arc and anode subdomains in the computational domain. However, the electrodes are assumed to remain solid and to the cathodes in the computational domain, in order to reliably account for the magnetic pinching induced by the tip of the cathode. On the other hand, the considered pointed cathode configurations are ideal geometries that can not exist in the real world. They lead to a geometrical singularity at the tip, that results in spurious computations in its vicinity, which overestimate the physical quantities. However, they enable to capture the trends up to the limiting ideal case. Finally, the general agreement between our results and the discussed experimental observations reinforces the validity of our numerical results, and confirms the need to incorporate the cathode into the calculation domain.

Introduction

In this chapter, the study and analysis of a fully coupled (cathode-arc-pool coupling) 2D axisymmetric TIG Spot configuration is presented. Here the interaction between the arc-plasma and the weld pool is of interest. To this end, the thermohydraulics of a 316L weld pool is considered and the behaviour of the TIG system is studied. To ensure that both the numerical and physical consequences of the arc-pool coupling problem are understood, both a numerical and physical discussion are presented. Effectively, the discussion begins, in section 5.6.1, with a comparison of the performance of the Dirichlet-Neumann and the quasi-monolithic algorithms implemented in the fluid model. The chapter follows with a brief discussion of the weld pool geometry and its behaviour as a function of inlet current, interface hypothesis and the chosen viscosity of liquid 316L steel, in subsections 5.6.2.1, 5.6.2.2 and 5.6.2.3, respectively. The analyses performed in this chapter serve to identify and study, at low computational costs, the dominant thermohydraulic effects at play, and the characteristics of the arc-pool system.

Geometric configuration

The geometric configuration used in this section is briefly discussed and presented. The schematic of the geometry is presented in figure 5.1. The study performed in this section is based on a TIG spot, axi-

Boundary conditions

The boundary conditions presented in table 5.1 are imposed onto the geometric configuration that is schematised in figure 5.1. A plasma inlet boundary condition is added at segment [FJ] which corresponds to a Q = 16 L•min -1 flow rate. Furthermore, the velocity symmetry condition is extended to encompass the anode domain. This serves as a symmetry condition for the pool domain as it forms during a calculation.

The imposed temperatures along the anode boundaries are chosen so as to maintain a cool workpiece and a moderately sized pool. As for the API at [BC], it no longer is considered as a boundary condition for the momentum equations, rather, an interface condition between the arc and the pool. The details of interface treatment are found in chapter 3.

Material properties

In this section, a solid thoriated tungsten cathode and argon arc are also used, however, unlike chapter 4, the anode is no longer considered to be a solid copper piece. With the focus of this section being the arc-pool dynamics of the system, the incorporation of a liquid melt pool material is paramount. To this end, the thermophysical properties of a 316L steel is chosen. The specific enthalpy, specific heat capacity, thermal conductivity, electrical resistivity and the mass density of the steel are extracted from Pichler et al. [START_REF] Pichler | Measurements of thermophysical properties of solid and liquid nist srm 316l stainless steel[END_REF]. The material properties are presented in appendix A. As for the viscosity of liquid steel, the values generally found in the literature on the thermo-hydraulic modelling of welding are higher than those reported in experimental studies. This is encountered in the works published by a number of authors, see [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Chien Nguyen | Sensitivity to welding positions and parameters in GTA welding with a 3D multiphysics numerical model[END_REF][START_REF] Traidia | Multiphysics modelling and numerical simulation of GTA weld pools[END_REF] for the numerical studies, and [START_REF] Choong | Thermophysical properties of stainless steels[END_REF][START_REF] Bobkov | Thermophysical properties of materials for nuclear engineering: a tutorial and collection of data[END_REF] for the experimental reports.

This section will explore some differences in the arc-pool thermal and dynamic responses as a function of different liquid steel viscosities. Thus, two dynamic viscosities are utilized in this section, with one extracted from the numerical work of Nguyen [START_REF] Chien Nguyen | Sensitivity to welding positions and parameters in GTA welding with a 3D multiphysics numerical model[END_REF], and the other from an IAEA nuclear materials data report [START_REF] Bobkov | Thermophysical properties of materials for nuclear engineering: a tutorial and collection of data[END_REF]. The viscosities used are listed below:

∀ T > T liq ; µ = 0.03 Pa • s (5.1)
Viscosity used by Nguyen [START_REF] Chien Nguyen | Sensitivity to welding positions and parameters in GTA welding with a 3D multiphysics numerical model[END_REF]

µ =                        0

Discretisation and meshing

This subsection discusses the chosen discretisations used for the different primal variables, the different meshes used and the spatial convergence of some key variables that the arc-pool model solves for. The discussed mesh dependent variables of interest are the pool temperature, velocity and pool contour geom-

Electric ϕ = 0 [HI] -σ∇ϕ • n = 0 [EH] ∪ [IC] ∪ [CD] ∪ [DF] -σ∇ϕ • n = j imp [EF] Heat transfer T = 750 K [HI] T = 500 K [EF] ∪ [FJ] ∪ [CI] ∪ [KD] -λ∇T • n = 0 [EH] ∪ [DC] ∪ [JK] Magnetic B • θ = 0 [EH] Momentum u = (0, 0) [JK] ∪ [FG] ∪ [GZ] ∪ [ZA] u • r = 0 [AH] ∪ [KD] u = (0, -1.4) m•s -1 [FJ] (σ • n) • r = 0 [DC] (σ • n) • ẑ = 0 [KD] ∪ [AH]
Table 5.1: Boundary conditions imposed along geometry contour seen in figure 4.1.

etry. Once an appropriate mesh is identified, it is used for similar configurations throughout the rest of this section.

Discretisation

With the velocity and interface deformation variables being the most mesh sensitive variables, they require more calculation nodes in a given mesh when compared to the other variables in the physical model.

Therefore, in the total domain, the velocity and interface deformation variables (u, v, h z ) are discretised using (Q2, P2) elements, while the voltage, temperature, magnetic field, pressure variables and the temperature depenendent physical parameters (σ, λ etc...) are discretised using (Q1, P1) elements. Assigning (Q2, P2) discretisations to the most mesh sensitive variables saves calculation time per global iteration when running simulations while ensuring sufficient spatial resolution of the solutions. The choices of element discretisation per physical variable and coefficient are presented and summarized in table 5.2.

The different meshes

The global mesh presented in figure 5.2 is made up of the cathode domain in red, the arc plasma domain in green and the anode domain in blue. The anode region, being the domain where the pool is formed, is meshed with a higher concentration of structured elements in the expected pool zone. The pool domain being of interest in this section, mesh refinement is concentrated in the refinement of the rectangular and structured mesh zone of the anode. Referring to figure 5.2, the structured mesh zone in the anode is succesively refined, until spatial convergence, using the following mesh densities:

ρ M = 8 3 , 4 3 , 2 3 , ... × 10 -4 m Interpolation Elements (ϕ, T , h, σ * , λ, etc...) → Q1, P1 (u, v, h z ) → Q2, P2

Spatial convergence of important variables

The spatial convergence of the different primal variables generally depends on multiple simulation settings, and in general, higher intensity configurations require finer meshes. To account for the mesh dependence of the simulations run for this section, a brief spatial convergence study is presented. The spatial convergence study is based on the most intense TIG Spot configuration (I = 150 A) studied in this section, chosen because of the sharp gradients in its results. The configuration has the most mesh sensitive solutions and so it ensures the sufficient spatial resolution of all other studied TIG Spot configurations. It is noteworthy to state that the arc variables are assumed to be sufficiently resolved with the chosen meshes because of the meshing choices taken in chapter 4. Furthermore, because this section is interested in studying the influence of different viscosities onto the arc-pool system, care is taken when treating the weld pools and claiming the spatial convergence of the simulations.

It is well known in computational fluid dynamics that the spatial convergence of a velocity field on a given mesh can depend on the Reynolds number of the flow. This can be understood by the increase in the gradients of the velocity fields as the Reynolds number of a flow increases; thus, an increase in mesh elements becomes necessary for the capture of the details of that flow [START_REF] Kuzmin | Finite element methods for computational fluid dynamics: a practical guide[END_REF]. This means that if for the same TIG Spot configuration the viscosity from equation (5.1) is chosen, the spatial convergence of the velocity fields occurs at a coarser mesh than if the viscosity from equation (5.2) is chosen. Moreover, with the temperature field in the pool being dependent on the pool velocity, the spatial convergence of the temperature field and the pool geometry is dependent on the viscosity value chosen. Therefore, the spatial convergence of the variables of the same TIG Spot configuration is studied but for the two different viscosities identified in subsection 5.4.

Having established the influence of the viscosity of the melt pool onto the spatial convergence of a simulation configuration, the following cases are studied and the resolution of their solutions analysed.

The cases are identical except for the different liquid steel viscosities studied.

Case 1:

I arc = 150 A α = 15 o r int = 0.
I arc = 150 A α = 15 o r int = 0.3 mm c S = 10 ppm µ → equation (5.2) Interface hypothesis → deformable
The analysis begins with case 1 which uses the higher viscosity value as was taken from the work of Nguyen [START_REF] Chien Nguyen | Sensitivity to welding positions and parameters in GTA welding with a 3D multiphysics numerical model[END_REF]. Three meshes are generated and used to simulate the configuration of case 1 above. The meshes are numbered from the coarsest to the finest, going from 1 to 3. The coarsest and finest meshes, mesh 1 and 3 respectively, are schematised in figure 5.2. The mesh densities of the structured zone in the anode go from ρ M1 = 8 3 × 10 -4 m to ρ M3 = 2 3 × 10 -4 m. From the simulation results of the different meshes, the different velocity and temperature fields are compared. The monitored variables of interest are presented in figure 5.3 and are the pool geometry, temperature at the API, and the pool velocity at both the API (magnitude) and the symmetry axis (z-component). Similar to section 4.5.3, the variables monitored for their spatial convergence are analysed with the relative difference equation. Comparing the results, the relative difference from equation (4.4) is used, and it drops below < 10% for the monitored variables. Thus, the results calculated on mesh 3 are considered spatially converged. Analogously, in looking for the spatial convergence of case 2, which uses the lower viscosity values as they were extracted from the IAEA report [START_REF] Bobkov | Thermophysical properties of materials for nuclear engineering: a tutorial and collection of data[END_REF], three meshes are generated and used for simulation. The meshes are numbered from coarsest (mesh 1) to the finest (mesh 3). To achieve spatial convergence of the solutions for case 2 is slightly more difficult in that finer meshes are required. In fact, the coarsest mesh used for case 2 is the finest of that used in case 1. The mesh density of the structured zone in the anode (see figure 5.2) for the coarse mesh 1 has a mesh density of ρ M1 = 2 3 × 10 -4 m, and the meshes are refined up to ρ M3 = 1 6 × 10 -4 m for mesh 3. Using the same variables of interest to monitor the spatial convergence of the simulations, the pool geometries, temperatures at the API, and the pool velocities at both the API (magnitude) and the symmetry axis (z-component) are plotted in figure 5.4. Comparing the results, the relative difference from equation (4.4) is used, and it drops below < 10% for the monitored variables. Thus, the results calculated on mesh 3 are considered spatially converged.

Discussion and results

This subsection discusses two topics. First, a brief discussion of algorithm performance is presented. The studied algorithms are concerned with the coupling techniques of the momentum equations between the arc and pool domains. Second, a brief discussion of the sensitivity of the pool domain to the electric intensity, choice of the interface hypothesis and the liquid steel viscosity values is presented.

Comparing algorithm performance

Here, a brief numerical experiment is set up, executed and discussed in an attempt to choose the appropriate coupling algorithm for the resolution of the momentum and mass equations in the unified model of this work. The numerical experiment and the following discussion are based on numerical tests comparing the partitioned Dirichlet-Neumman to the quasi-monolithic algorithms implemented in this work. To be able to compare the performance and results of the two implemented algorithms, the simulations need to be set up in an analogue manner to one another. To this end, eight cases are set up, run and then analysed. Four setups are identifed and are based on three TIG Spot configurations where for one of the configurations the relaxation factor is varied. The four identified setups are executed using the two implemented coupling algorithms. To be consistent in the comparisons between the different cases, all tests are studied on the same mesh. The mesh is based on one anode mesh density, with a value at ρ M = 1 3 ×10 -4 m. Furthermore,

c S = 10 ppm µ → equation (5.2) Interface hypothesis → deformable                  ) I = 100, 125, 150 A (5.3)
The three configurations listed above are run on both algorithms for the same relaxation factor applied to the momentum and mass equations, at a value of ω N S = 0. solid. The other cases are analogously initialised.

Convergence criteria and linear system solvers

To ensure that the simulated cases are sufficiently converged when comparing their results and performance together, the convergence criteria used in the algorithms are lowered. The algorithm convergence criteria specific to the performance study of this section are summarized in the table below. Furthermore, since

δϕ min δT min δu min δv min δh z min 2D 1 × 10 -4 1 × 10 -6 1 × 10 -6 1 × 10 -6 1 × 10 -6
Table 5.4: Convergence criteria used only in this chapter.

this section is only concerned with the performance of the algorithms used in the fluid dynamics module, only its linear solver is of concern. To ensure that the resolution of linear systems that have either well or ill-conditioned matrices is performed in a comparable manner, the direct solver programmed in Cast εM is used. The solver is based on an LDL t method for symmetric matrices and the Crout method for nonsymmetric matrices. The performance of these methods depend on the linear system being treated and on the machine executing the program, but they are of the order of:

speed ∼ O(N 4/3 ) (5.4)
where N is the dimension of the linear system [START_REF]Cast3m Web site[END_REF]. Finally, it is important to note that the computer used to run the simulations has a quad-core Intel i7-7700 CPU with a baseline clockspeed at 3.60 GHz and 8.192, 0.256 and 0.032 Mb for the L3, L2 and L1 cache memories. The computer also houses 4 slots of 8

Gb DDR4 RAM chips that operate at a clockspeed of 3200 MHz.

Quantifying the performance

In order to quantify the differences in performance between the two algorithms, both the total calculation time and the total number of global iterations of a given case are used. Furthermore, to verify that both algorithms converge to solutions that are close to each other, key variables are monitored and used for comparison. Also, the increments of the different primal variables are also presented, as functions of iteration number, in algorithm convergence plots. Thus, the sub-cases (a and b) of each case from table 5.3 are compared to one another. The analysis begins with a presentation of the differences in both the convergence plots and their solutions. Starting with case 1 the following key quantities and variables are presented: From the presented results, it is apparent that for case 1, the quasi-monolithic algorithm performs slightly faster and requires a smaller number of iterations to achieve convergence. From figures 5.5 c) and d), it is noted that the two algorithms give simulation results that are almost identical. From the presented results, it is apparent that for case 2, the quasi-monolithic algorithm performs slower and requires a larger number of iterations to achieve convergence. From figures 5.6 c) and d), it is noted that the two algorithms give simulation results that are almost identical. From the presented results, it is apparent that for case 3, the quasi-monolithic algorithm also performs slower and requires a larger number of iterations to achieve convergence. From figures 5.7 c) and d), it is noted that the two algorithms give pool velocity profiles that are slightly different but the difference can be considered negligible because it insignificantly influences the pool geometry. From the presented results, it is apparent that for case 4, the Dirichlet-Neumann algorithm does not achieve convergence. The calculation is stopped at global iteration (k = 2000) as it is evident that no stationary solution was found. As for the quasi-monolithic algorithm, it converges at iteration (k = 948). It is important to note that in referring to figures 5.5, 5.6 and 5.7, that the algorithms are mathematically consistent as they both tend towards the same solutions at convergence. Furthermore, the different performances of the cases presented reveal that a relaxed Dirichlet-Neumann algorithm converges quicker than the quasi-monolithic algorithm. The total calculation times and iteration number until convergence (or until the calculation is stopped) for all the studied cases are summarised in table 5.5. While an in-depth algebraic analysis of the results on algorithm performance presented above is out of the scope of this thesis, the following remarks are discussed based on figures 5.5, 5.6, 5.7, 5.8 and table 5.5:

Remark 1 : It is apparent that between case 1 and case 3, the total run time and total number of iterations until convergence increases. The increase in the electric current that passes through the system makes the unified TIG model more difficult to solve, independently of the algorithms. This is because as the arc intensity rises with a rise in inlet electric current, the stationary solution to the configuration gets farther away from the field at initialisation 2 .

Remark 2 : For the 2D configurations treated in this section, the direct solver shows negligible discrepancy in execution time per iteration between the relaxed linear systems of the Dirichlet-Neumann and the quasi-monolithic algorithm.

Remark 3 : It is apparent that the relaxed Dirichlet-Neumann algorithm achieves convergence faster than the relaxed quasi-monolithic algorithm between cases 1 and 3 because it requires a smaller number of iterations to convergence. However, in referring to case 4.b), if the quasi-monolithic algorithm is executed without relaxation, the simulation requires a smaller number of iterations before convergence. However, the lack of relaxation makes the direct solver require more time to solve the linear algebraic system at every global iteration of the simulation.

Remark 4 : Comparing case 4.b) to 2.a), the un-relaxed quasi-monolithic algorithm requires a smaller number of total iterations but a longer total calculation time until convergence than the relaxed Dirichlet-Neumann algorithm. Quantifying the remark, the relaxed Dirichlet-Neumann algorithm requires 22% more iterations but 29% less calculation time until convergence. This is because solving the un-relaxed linear system using the direct solver requires about 1.7 times more time than the relaxed linear system.

Remark 5 : Comparing case 4.b) to case 2.b), the relaxed case, although requires a larger number of iterations for convergence, converges faster than the un-relaxed case. Quantifying the discrepancy, case 2.b requires 46% more total iterations for convergence, but 13% less total calculation time until convergence. This is because solving the un-relaxed linear system using the direct solver requires about 1.7 times more time than the relaxed linear system.

Remark 6 : Referring to case 4.a), it is apparent that the lack of relaxation of the Dirichlet-Neumann algorithm makes that it does not converge to the required convergence criteria, rather it displays an oscillatory solution. This is because the force balance at the interface of coupling is not satisfied with an average of 2.5 internal iterations per global iteration. The algorithm ends up over and under shooting the force balance at the coupling interface, thus preventing the convergence of the algorithm.

The force imbalance of the Dirichlet-Neumann algorithm is hypothesized to be related to the ratio of the numerical Reynolds numbers in the two coupled domains. A formal study of this hypothesis is reserved for a future work.

Following the above remarks, it can be argued that even though monolithic schemes are known to have no instabilities associated with the interface treatment [START_REF] Michael | Stability analysis of numerical interface conditions in fluid-structure thermal analysis[END_REF] Furthermore, in going towards a 3D approach, the quasi-monolithic algorithm is expected to perform significantly slower than the Dirichlet-Neumann algorithm. This is because a direct solver does not scale linearly with the number of degrees of freedom in a linear system3 . Finally, the use of the Dirichlet-Neumann algorithm must be performed with sufficient relaxation and internal iterations so as to achieve a satisfactory force balance at the interface. Thus, a maximum of 3 internal iterations is allowed in the algorithm and a relaxation factor of ω N S = 0.3 is used for the rest of this work.

Pool sensitivity to physical parameters and interface hypotheses

The weld pool is generally influenced by a multitude of parameters, including the arc electric intensity and the viscosity of the liquid metal [START_REF] Mills | Factors affecting variable weld penetration[END_REF]. In this section the influence of the arc intensity, the liquid steel viscosity and the interface hypothesis are discussed. These parameters are of interest because they engender the influence of both arc and pool parameters onto the response of the weld pool. The analysis begins by discussing the influence of the imposed electric current onto the weld pool, then follows with a discussion of the influence of the interface hypothesis onto the pool and finally the influence of the viscosity of steel onto the pool. The significant rise in the velocity of the weld pool with a rise in the inlet electric current increases the rate of heat advected away from the pool center. These effects influence the pool size and its depth to width ratio. The weld pools calculated for the different inlet electric currents are presented in figure 5.13. Furthermore, the pool streamlines for the cases at I = 100 A and I = 150 A are presented in figures 5.12.a) and b), respectively. The stream function in a cylindrical coordinate system is defined as:

ψ = 2π L r(u • n)dL
where n is the normal vector to an arbitrary line element dL. Reminding the reader that the stream function is defined as the volume flow rate that passes through the arbitrary line element dL. This measure allows to visualize the streamline of a flow, which facilitates analysis. So, in referring to the pool geometries in

Pool width:

In this work, the pool width is defined as the pool maximum radius. As is apparent from figure 5.13, an increase in the inlet electric current increases the width of the weld pool. In fact, the width of the pool increases from 5.95 mm to 9.14 mm between the I = 100 A and I = 150 A pool; which is a 1.5-fold increase. Analysing the forces responible for this growth, I start by noting that according to the work of Nguyen, neither the bouyancy force nor the Lorentz force in the pool are the dominant forcing terms [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF].

Furthermore, for pools at low concentrations of sulfur, the Marangoni force4 although it forces the pool at the same order of magnitude as the shearing of the arc, is not as variable w.r.t inlet current. This is visualized in figure 5.14, where the Marangoni forcing increases negligibly when compared to the rise in . It becomes evident that as the arc shear stresses increase, the pool surface velocities consequently increase. The increase in the surface velocities, along with the increase in transported heat from the arc, increases the advected heat towards the periphery of the pool. Thus, a rise in electric current increases both the shearing and the heat flux transferred to the pool surface which increases the advected heat outward and away from the center, which widens the pool.

Pool depth:

As is apparent from figure 5.13 the depth of the pool increases with a rise in the inlet electric current. In fact, the depth of the pool increases from 1.76 mm to 3.22 mm, which is a 1.8-fold increase. Although the depth of a weld pool increases in a similar manner to the pool width, the pool is wider than it is This is due to the weak downward velocities, which makes the pool penetration in these cases dominated by simple diffusion, which is apparent from the quasi-constant descent of the temperature profile along the pool symmetry axis. However, the deepening of the pool near the center for the I = 150 A case is dictated by a secondary circulation zone, mapped in red/yellow in figure 5.12. This secondary circulation zone exists for all cases but is significantly stronger for the I = 150 A case. Moreover, the secondary circulation is identified to be a Lorentz eddy because if the Lorentz force is removed from the weld pool model, the eddy is not calculated. This effect is presented in figure 5. [START_REF]Cast3m Web site[END_REF] where the I = 150 A configuration Consequently, the pool remains wider than it is deeper for all cases because the secondary circulation zone is significantly less intense than the primary circulation. However, the more intense the inlet current, the higher the intensity of the secondary circulation; and with the Lorentz eddies flowing in a counter-clockwise manner, they advect the hot liquid downward. Quantifying the pool geometries table 5.6 presents the widths, depths and the depth to width ratio of the three weld pools. The calculated D:W ratios show that the deepening effect of the secondary circulation eddy is 1.5-fold higher between the I = 125 A case and the I = 150 A case. This ratio increase is hypothesised to be due to a strengthening of the interaction of the primary circulation with the Lorentz eddy. Effectively, the primary circulation and the Lorentz eddy intensify with a rise in the inlet electric current, which strengthens their interaction.

Once their interaction is strong enough the downward advection of hot liquid is significantly increased (as seen in the I = 150 A case), augmenting heat transport locally. This consequently contributes to a stronger deepening of the pool, with the deepening concentrated near the center and below the red circulation zone from figure 5.12.b).

For purposes of comparison, the experimental studies reported by Mills [START_REF] Mills | Factors affecting variable weld penetration[END_REF] show a general agreement with the trends discussed in this subsection. Indeed, a rise in the inlet electric current leads to an increase in the experimental dimensions of the weld pool.

Influence of the interface hypothesis

Here, I briefly discuss the influence of the interface hypothesis onto the pool system. Although the Weber number of the arc-pool system is generally small (see chapter 1, table 1.1), the choice of the interface hypothesis can significantly affect the system. This is due to the strongly coupled thermo-hydraulic nature of the system. The discussion is based on the following two configurations: By maintaining all the simulation parameters constant except for the interface hypothesis, the effect of the weak deformation of the interface onto the pool system can be discussed. Figure 5.19 presents the pool geometries of both case 1 and 2. As the pool contours suggest, allowing for a deformable interface in the arc-pool system increases both the width and depth of the pool. However, the pool width increase is negligible (≈ 2%) as compared to the pool depth increase (≈ 17%). The significant increase in depth is primarily due to the distortion in the shape of the melt pool. The deformation at the interface modifies the location and intensity of the secondary circulation below the interface and thus also the heat advected downard. Figure 5.17 shows the streamfunction field for both case 1 and 2 in the vicinity of the secondary circulations. The eddy for case 2 is seen to be larger and extends further downward as compared to that seen for case 1. This can be explained using the following two arguments. Due to the incompressibility of the melt pool the depression of the interface shifts the liquid steel that is directly underneath it ≈ 0.24 mm deeper into the anode. This downward displacement of the weld pool modifies the geometry of the basin the eddy is placed as seen in figure 5.17. Furthermore, eddies are known to be sensitive to the boundary conditions of their basin [START_REF] Moffatt | Viscous and resistive eddies near a sharp corner[END_REF] and so the modification in the pool geometry appears to intensify the eddy computed in case 2. This intensification increases the advected heat downward into the anode and so extends the penetration of the pool an additional ≈ 0.22 mm.

Influence of the liquid steel viscosity

Here, I briefly discuss the influence of the choice of the viscosity of the liquid steel onto the pool system.

With weld pools being mainly surface driven fluid systems, a variation in the dynamic viscosity will translate to the surface velocity, which then translates to the thermal and the secondary circulations that are dependent on the surface flow. The discussion is based on the following two configurations:

Case 1: By maintaining all the simulation parameters constant except for the choice of the dynamic viscosity data set used, the effect of viscosity onto the pool system can be discussed. Therefore, by maintaining all arc parameters constant it is negligibly affected by the change in the viscosity of the pool. This is seen in figure 5.18 a) and b), where neither the arc velocity field nor the temperature field is significantly affected.

I arc = 150 A α = 15 o
Consequently, both case 1 and 2 have arcs that impart similar dynamic stresses and transfer similar heat fluxes to the pool surface. However, the computed pool geometries are significantly different in both their width and depth. This is presented in figure 5. [START_REF] Cuvelier | Some numerical methods for the computation of capillary free boundaries governed by the navier-stokes equations[END_REF] where the pool width varies 5% between cases and the depth is computed to be 1.8mm for case 1 and 3.2 mm for case 2, which is a significant 1.7-fold difference.

With the arc remaining almost invariable between the cases, the computed differences in the pool geometry are due to the thermo-hydraulics of the pool. The identified effects driving the sensitivity of the pool to the dynamic viscosity are threefold, and are based on figure 5.20:

Effect 1 : The pool at a higher viscosity is stiffer than that at a lower viscosity and so for the same arc shearing results in lower surface velocities. The surface velocity norms for cases 1 and 2 are plotted in figure 5.20.a) where the velocity profile rises 4-fold between case 1 and 2. The higher surface velocities slightly increase the width of the pool.

Effect 2 : The higher velocities in the low viscosity pool allow for a more significant advection of heat around the pool, that which allows for pools with higher internal energies at the equilibrium stationary solution. The more easily internal energy is distributed in the pool, the more internal energy can be stored in the pool at its thermal equilibrium state. This effect is seen in figure 5. 20.b) where the center of the pool is 1.15-fold hotter for the less viscous case.

Effect 3 : Due to the weaker surface velocities for the more viscous pool, its consequent primary circulation zone is also weaker. The higher viscosity also stiffens the pool which dampens the effect the Lorentz force has on the pool. Furthermore, the lower pool temperatures reduce the overall volume of the pool, that which does not allow for a secondary circulation to form. This is seen in figure 5.20.c) where for the more viscous pool, there is no downward flow along the symmetry axis, which is associated to the secondary circulation. The lack of the downward advection of heat means that the penetration of the high viscosity pool is diffusion dominant, which can be deduced from figure 5.20.d). Thus, with lower pool temperatures, and no secondary circulation, no heat is advected downward which explains the 1.7-fold difference in pool depth between case 1 and 2.

The sensitivity of the pool to the choice of dynamic viscosity requires that care is taken when setting choosing the model parameters. Although lower pool viscosity values requires that the model be solved on a finer mesh (see section 5.5.3), which makes for more costly computations, the consequent physical

Nota bene:

It is important to note that the thermal boundary conditions can strongly influence the behaviour of both the arc and pool. However, they are expected to only influence the magnitude of the variables and not the captured trends. An in-depth analysis of this influence is out of the scope of this thesis.

Summary

Studying the influence of the coupling algorithms onto the arc-pool model ensured that the more performant method is identified and used in the rest of this thesis. Indeed, the two studied algorithms were found to converge to very close solutions and thus were considered comparable. Furthermore, in comparing time to convergence and total iteration number, the Dirichlet-Neumann coupling algorithm was found to achieve numerical convergence faster than the quasi-monolithic algorithm. Additionally, the Dirichlet-Neumann algorithm was argued to be more appropriate for use in 3D configurations. It is interesting to note that based on remarks #4 and #5 in subsection 5.6.1.2, there appears to be an optimal relaxation factor that would reduce both calculation time and iteration number to a minimum.

Studying the influence of the inlet electric current, the interface hypothesis and the liquid steel viscosity value, at low computational costs in a 2D axisymmetric configuration, served to identify the dominant thermohydraulics at play in the arc-pool system. In fact, the three model parameters were found to significantly influence the calculated dimensions of the pool. Effectively, a rise in inlet current increases both the arc shear force at the API and augments the strength of the secondary circulation in the weld pool. These effects were shown to increase both the width and depth of the pool. Furthermore, although the Weber number of the arc-pool interaction is small, the deformation of the API considerably influences pool penetration by both displacing and morphing the secondary circulation of the weld pool. As for the influence of the pool viscosity onto the arc-pool system, the main effects are the following. Using an augmented pool viscosity lowers the maximum velocities and temperatures calculated in the pool which lowers the widening and penetration of the pool. Moreover, the augmented rigidity of the pool strongly counteracts the Lorentz force and obstructs the formation of a secondary circulation. This further reduces the penetration effect of calculated weld pools.

Introduction

A presentation of two different 3D welding configurations is analysed and discussed in this chapter. The chapter begins with section 6.1, where a presentation of a verification study performed on a TIG spot configuration modelled in 3D Cartesian coordinates is discussed. The purpose of the study is to verify that the extension of the 2D axisymmetric model to a 3D one was coherently performed. Next, in section 6.2, the verified 3D TIG welding model is extended to account for weld displacement effects. The extended model is then used to set up a preliminary validation study performed against an experimental case chosen from the literature. Some results of the simulated case are presented and the thermohydraulics of the welding configuration briefly discussed. Finally, both similarities and discrepancies between the numerical test and the experimental case are analysed and discussed; moreover, the limitations of the validation study are presented.

A TIG spot verification study

In this section a verification study of the 3D model developed in this work is presented. This verification study is based on the assumption that rotational symmetry exists for steady state spot configurations;

thus, the verification study is performed on a spot configuration where a 2D axisymmetric configuration is compared to its 3D Cartesian analogue. The 2D axisymmetric spot configuration is set up, executed and its simulation results processed so as to serve as a reference solution to compare the simulation results of the 3D analogue to. The analogous 3D spot model is set up using the same welding parameters as the reference 2D model, but is built on a full 3D Cartesian mesh. Furthermore, the 3D spot configuration is solved on a coarse and a fine mesh, which allows for a brief discussion of spatial convergence. To begin the discussion, first the geometric configuration of the study is presented.

Geometric configuration

Being constructed on analogous geometries, the 2D axisymmetric and 3D Cartesian spot configurations are modelled using one schematic. The schematic is presented in figure 6.1, where the faded dash-dot symmetry axis represents the axis of rotation. The schematic represents a z -x section of the geometry centered at y = 0. Furthermore, the 2D axisymmetric model is built on the half geometry of the schematic (see figure 5.1 in chapter 5), while the 3D Cartesian model is built on the full rotation of the schematic.

Boundary conditions and material properties

The boundary conditions imposed onto the reference 2D axisymmetric model are the same as those presented in table 5.1 from chapter 5; while the boundary conditions imposed onto the 3D spot model are presented in table 6.1. The 3D spot boundary conditions are presented as they apply to the schematic in 6.1 and 6.2 of this chapter.

Discretisation

To reduce computational costs when simulating the 3D model all the solved variables are discretised tri-linearly (u, P, v, p, h z , h, B, ϕ) while the temperature dependent physical parameters are constant by element. Furthermore, the volume of the geometry Ω is discretised using a mix of hexahedral and prism elements, while the surface and interface subdomains ∂Ω, Γ are discretised using a mix of quadrilateral and triangular elements. The choices of element discretization per physical variable and coefficient are presented and summarized in table 6.2. The tri-linear Q1-Q1 and P1-P1 discretisation of the velocity and pressure variables in this chapter does not respect the Ladyzhenskaya-Babuska-Brezzi condition. Thus, to ensure the stability of the coupled velocity-pressure scheme, the system is stabilised with polynomial pressure projection terms [START_REF] Clark R Dohrmann | A stabilized finite element method for the stokes problem based on polynomial pressure projections[END_REF].

Meshing

Two 3D meshes are generated on a Cartesian coordinate system for this verification study and they serve to simulate the 3D spot model in a manner analogue to the 2D axisymmetric spot model. The meshes are constructed by rotating a 2D axisymmetric mesh (see figure 5.2 in chapter 5) 360 o around its symmetry axis. The rotation of the 2D plane is discretised with N rota steps, which generates the discretised third spatial dimension. The two generated meshes, mesh 1 and mesh 2 have N rota = 15 and N rota = 30 steps, respectively. Figure 6.2 schematises the generation of the two meshes used in this section. Figure 6.2: Mesh generation of 3D spot model for two rotational discretisations.

Discussion and results

To solve the stationary TIG spot model on a 3D mesh, spot configurations with steady-state equilibrium solutions should be used, otherwise the simulation would risk diverging. Thus, to simulate a spot configuration that easily converges and allows for the verification of the 3D model, the augmented viscosity value adapted from Nguyen [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF] is used. The main weld parameters of the chosen spot configuration are listed below:

I arc = 150 A c S = 0 ppm µ = 0.03Pa • s Interface hypothesis → deformable.
The above configuration is simulated on both mesh 1 and mesh 2, and the results of the converged simulations are discussed in two parts. First, the results of the 3D spot model on the finer mesh are presented, analysed and discussed, with a focus on the rotational symmetry of the solution. The second part of this verification study compares the results of both the coarse and fine meshes to the simulation results of their analogous 2D axisymmetric model (run using the same spot configuration). It is important to note that the verification of the numerical fluid dynamics and magnetostatic schemes were performed in the work of Brochard [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF]; moreover, the extensive verification of the used mathematical operators was performed in Gounand's work [START_REF] Gounand | Modèles multiphysiques arc-bain en soudage à l'arc TIG[END_REF].

Comparison of the 2D to the 3D model

To compare the 3D model to the reference 2D axisymmetric model, the temperatures and velocities of the arc-plasma and the weld pool, and the weld pool geometries are compared. The simulation solutions are generated on mesh 1 and mesh 2 of the 3D spot model, and the reference solution on the analogous 2D axisymmetric mesh (see figure 6.2). The temperature and velocity variables across the arc symmetry axis, the API and the pool symmetry axis are plotted in figure 6.6. As observed in the plotted results, the 3D spot model reproduces the same results and trends calculated by the 2D axisymmetric model. In fact, the temperature fields calculated with the 3D model reproduce, to a good level of precision, the field calculated with the 2D model. The velocity fields calculated with the 3D model also reproduce the trends calculated by the 2D model; however, the precision is slightly lower. Furthermore, the three calculated pool contours show good agreement as is seen in figure 6.7. Although the velocity field calculated on mesh 2 is in better agreement with the reference 2D field, mesh 1 is considered satisfactory because its relative difference w.r.t the reference solution remains < 10% (see (4.4)). In fact, the small relative differences in the plotted profiles can be attributed to numerical integration errors and are expected to reduce towards 0% as the 3D mesh is further refined, and in the limit that the algorithm convergence criteria are reduced to 0.

Thus, both the rotational symmetry of the 3D spot model, and its good agreement with the 2D axisymmetric spot model imply that the extension of the 2D axisymmetric model to a 3D Cartesian model is mathematically consistent. Moreoever, the algorithms defined in chapter 3 are functional and conservative in both 2D and 3D configurations.

TIG welding with displacement effects

In this section, the 3D model is extended to account for displacement effects that occur in a typical welding operation. The displacement effects are accounted for fixing the cathode w.r.t the laboratory frame of reference and assigning a weld displacement velocity to the workpiece. To test the extended model, a welding configuration is chosen from Koudadje's thesis [START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF] and is simulated. The thermohydraulic fields of the simulated welding configuration are discussed and analysed, and the weld pool geometry is compared to the experimental results from Koudadje [START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF]. The limitations of the preliminary validation study are also discussed.

Geometry configuration

To both account for the inlet/outlet heat flux associated with the weld displacement effects and to reduce computational costs, a new geometric configuration is designed. The geometry is set up with both anode inlet and outlet surfaces that allow for appropriate thermal boundary conditions to be applied. Furthermore, a z -x symmetry plane is assigned, which reduces the computational domain by half. The geometric configuration is schematised in figure 6.8. The workpiece thickness and arc height are chosen to correspond to the chosen experimental configuration in Koudadje's work [START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF]. However, the domain length and width are chosen to allow for sufficiently far inlet and exit surfaces. As seen in figure 6.8, the workpiece surfaces are aptly named ∂Ω anoin and ∂Ω anoout .

Boundary conditions, displacement effects and meshing

Here the boundary conditions used to close the 3D weld displacement model, the modelled displacement effects and the mesh used to simulate the welding configuration are presented. As for the material properties used in this section, they are the same as those presented in section 6.1.2.

Boundary conditions

The boundary conditions imposed onto the geometric model are similar to those of the previous sections, but are supplemented to account for the weld displacement effects. Furthermore, to capture displacement effects while using a steady state model requires that the model emulate welding on an infinitely long and wide but a finitely deep workpiece. To this end, a cold 300 K condition 1 is assigned to the inlet ∂Ω anoin surface, while the outlet ∂Ω anoout surface is left as a free -λ∇T • n = 0 W exit condition. Moreover, the workpiece bottom surface ∂Ω anobot is assigned a convection boundary condition with

h conv = 15 W•K -1 •m -2
and T inf = 300 K. All boundary conditions applied to the 3D weld displacement model are presented in table 6.3 as they apply to the geometry schematised in figure 6.8. As for the interfaces, the arc-pool interface is treated using the Dirichlet-Neumann algorithm (see section 3.6.2.1 in chapter 3), and the cathode-arc interface is considered to be a fixed, stationary surface where U wall = 0 m•s -1 .

1 Which emulates a workpiece that is at room temperature at infinity.

Meshing

One mesh is generated for the weld displacement model, and is meshed in three steps. First a 2D axisymmetric mesh (see figure 5.2 in chapter 5) is rotated 180 o and discretised using N rota = 30 steps. Second, a cylindrical region around the cathode and the central z -y plane at x = 0 are fixed. Lastly, the submesh in the right hand quadrant is compressed and the submesh in the left hand quadrant is extended to form the final mesh. The meshing procedure is schematised in figure 6.10. To ensure that the compressed and Figure 6.10: Mesh generation of 3D weld displacement model.

extended elements maintain their quality and orthogonality the mesh regularisation method adapted from

Huang [START_REF] Huang | Mathematical principles of anisotropic mesh adaptation[END_REF] is used.

Discussion and results

To test the capabilities and validity of the weld displacement model developed in this work, a numerical configuration analogous to a line weld case from the experiments studied in Koudadje [START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF] is set up. The welding configuration is adapted from Appendix F, weld line tests, experiment #2 from Koudadje [START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF] and is as follows:

Experimental parameters

)                                  I arc = 200 A H arc = 3 mm α = 20 o r int = not specified c S = 10 ppm H ano = 18 mm V displ = 30 cm • min -1 Workpiece material → 304L
where the cathode truncation radius is not reported. The numerical configuration that is set up with the above parameters except that the truncation radius is chosen at r int = 0.6 mm and the thermophysical properties of the workpiece material are for 316L steel2 . Furthermore, in order to ensure the simulation of a converging numerical configuration, without incurring excessive computational time 3 , the used liquid steel viscosity is based on the augmented viscosity value in equation (5.1). In fact, finding a converging solution to 3D configurations with the lower experimentally reported viscosity values (see equation (5.2))

takes a significantly longer amount of time and many more iterations. Moreover, only the results of a fixed interface simulation are presented because, while writing this section, the deformable interface simulation

had not yet converged.

Simulation results

Here I present the results of the simulation of the chosen welding configuration and briefly discuss the The temperature field of the global domain is calculated to be continuous across the arc-plasma and workpiece domains, as is represented by the black isocontours. The continuity of the temperature field implies that both the energy and momentum conservation laws are respected by the coupling algorithms implemented in this work (see sections 3.5 and 3.6). This indicates that the algorithms are robust and generally applicable to 3D welding configurations. Furthermore, the temperature field of the arc-plasma displays an asymmetry in its interaction with the workpiece and this indicates the importance of modelling the arc-plasma in 3D. Furthermore, the symmetry condition applied on the z -x plane can be verified by noting the orthogonality of the temperature isolines at the ∂Ω sym boundary. To this end, the temperature field at the arc-workpiece interface (Γ AP I ) is shown in figure 6.12. Effectively, the symmetry in the tem- Next, the arc-plasma flow, the hot regions of the arc and the workpiece, and the pool hydrodynamics are represented in figure 6.13 by the streamlines, temperature isosurfaces/contours and the vector fields, respectively. Indeed, the simulation reproduces the inlet argon flow and its acceleration and heating downstream the cathode which consequently generates the thermohydraulics of the weld pool, that in turn shapes the pool geometry. As is seen in figures 6.11 and 6.13, the displacement velocity of the workpiece stretches the temperature field of the arc-pool system and consequently breaks any rotational symmetries of the simulated system. The stretched pool geometry and its temperature field are represented in figure The incorporation of the cathode, arc-plasma and the anode/workpiece into the computational domain is important because it removes the need to experimentally calibrate equivalent heat and pressure sources and/or use 2D → 3D projection methods as was done by Nguyen [START_REF] Chien Nguyen | Sensitivity to welding positions and parameters in GTA welding with a 3D multiphysics numerical model[END_REF][START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF], Koudadje [START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF] and Traidia [START_REF] Traidia | Multiphysics modelling and numerical simulation of GTA weld pools[END_REF]. However, to ensure that the weld displacement model developed in this work is adapted to simulate experimental configurations, a validation study is necessary. Limited by time, a preliminary validation study is presented in the following section where the simulation results are compared to the experimental pool geometry as reported in Koudadje's thesis [START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF].

Comparison to experimental observables

To compare the calculated pool geometry to the experimental macrographic image reported by Koudadje [START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF], the appropriate calculated pool geometry projection must be considered. In fact, referring to figures 6.14 and 6.15, the pool geometry is not deepest where it is widest; furthermore, the pool floor exhibits multiple inflection points. Thus, calculating the equivalent macrographic slice from a stationary pool model requires that the maximum depth isocontour be calculated as a function of the x and y axes. The calculated maximum depth isocontour for the pool geometry is presented in figure 6.16 and is represented by the red points and contour on the mesh of the pool geometry. Using the equivalent numerical macrographic projection the results are compared to the experimental cross section in figure 6.17. Indeed, the calcu- Figure 6.17: Comparison of numerical maximum depth isocontour to the experimental macrographic cross section from Koudadje [START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF].

emission; d) the deformable interface hypothesis was shown in both subsection 5.6.2.2 and by Nguyen [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF] to influence the weld pool geometry; e) the thermal contact points of the experimental workpiece and the workbench, used in Koudadje's work [START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF], could influence the rate of cooling of the workpiece. Furthermore, it is possible that assigning a displacement velocity to the workpiece as opposed to displacing the welding torch (as done in Koudadje [START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF]) w.r.t the laboratory frame of reference, might influence the cooling rate of the arc and thus influence the pool geometry. In fact, Gonzalez [START_REF] Gonzalez | Numerical modelling of an electric arc and its interaction with the anode: part II. the three-dimensional model-influence of external forces on the arc column[END_REF] showed that the thermo-hydraulics of the arc is influenced by the cross-flow while displacing the torch; however, the cross-flow effect is only significant when the displacement velocity is high.

As a final note, the author encourages the development of an extensive validation study in the future, in order to thoroughly quantify the limitations and the extent of the validity of the developed model.

Summary

Spot verification study

The spot verification study ensured that the extension of the 2D axisymmetric model to a 3D Cartesian one was coherently performed. In fact, the calculated rotational symmetry of the 3D spot simulation indicates that the mathematical operators used to model the volume, surface and interface forces were properly extended from the 2D axisymmetric model. Moreover, the rotational symmetry of the 3D spot implies that the axisymmetric hypothesis in the 2D model is valid (for the spot regimes studied); thus, allowing for the use of the 2D spot model as a reference.

The comparison of the 3D simulation results to the reference 2D model ensured that extending the model did not modify the implemented physics in an unwanted manner. Furthermore, the strong agreement between the reference 2D axisymmetric and the analogous 3D Cartesian simulation results further indicates the mathematically coherent extension of the model. Finally, the agreement between the 2D and 3D models indicates that the coupling algorithms implemented in this work maintain their robustness and conservative nature in 3D configurations.

Welding with displacement effects

The weld displacement study incorporates the cathode, arc-plasma and the workpiece domains and presents the effects of a displacement velocity onto the welding system. The developed model is a potentially powerful tool that alleviates the need for the use of experimentally calibrated heat and pressure sources. As for the calculated asymmetry in the temperature field of the arc-plasma, it implies the importance of the use of a coupled arc-plasma and workpiece model when capturing 3D displacement effects. Furthermore, the continuity of the temperature field across the API indicates that the implemented coupling algorithms remain conservative when applied to weld displacement configurations.

The comparison of the numerical results to the analogous experimental case brings into context the validity of the model. The agreement between the calculated and measured pool depth implies that the model is capable of simulating experimentally measured trends. Furthermore, the reasons behind the weaker agreement between the calculated and measured widths are hypothesised. Although an exact agreement between the simulated pool geometry and the measured dimensions is not found, the fact that the simulation results are in same the order of magnitude as the experimental measurements implies the general validity of the model. Furthermore, the calculated pool temperature and velocity fields are similar to those encountered in the literature on the numerical modelling of weld pools [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF][START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF][START_REF] Traidia | Multiphysics modelling and numerical simulation of GTA weld pools[END_REF].

Conclusions and Perspectives

Conclusions

In this doctoral thesis, I focused on the coupling of previously decoupled numerical models for the magnetothermo-hydrodynamic simulation of TIG arc welding. This work presents the first fully coupled 3D cathodearc-pool welding model at CEA that accounts for both electrodynamic and thermohydraulic phenomena. A focus on the stability and robustness of the coupling algorithms and their applicability in both 2D and 3D models is privileged; moreover, the mathematical consistency of the developed and implemented numerical model is assured. The 3D model this thesis presents is sufficiently robust and developed to allow for the numerical simulation of welding with torch displacement effects.

Building on previous work, first, the electrodynamic, energy, mass and momentum laws were presented and the model hypotheses discussed. Care is taken to identify the dominant thermohydraulic contributions at the arc-pool interface; and the physical model is expressed in a form applicable to both 2D and 3D formulations. The physical model is then numerically implemented with a focus on multiphysics coupling strategies and the algorithms that handle the different thermo-hydraulic phenomena that dominate the arc-plasma and weld pool domains. Next, the sensitivity of the model to both physical parameters and to the implemented coupling algorithms is studied on 2D axisymmetric configurations. The analyses performed in 2D allow for a good trade off between acceptable physical representativity and affordable computational costs. After having explored the sensitivities of the welding model, a verification study ensures that no mathematical inconsistencies were introduced into the fully coupled 3D model. The verified 3D model is then used to study a welding configuration that accounts for weld displacement effects. This 3D weld displacement model is closer to industrial configurations. Thus, one of the few existing studies on 3D multiphysics modelling of TIG welding using a fully coupled cathode-arc-pool model was introduced.

In chapter 2, the physical model describing the dominant magneto-thermo-hydrodynamic phenomena present in TIG welding is discussed. The dominant electrodynamic, conjugate heat transfer, phase change and fluidic phenomena are mainly modelled based on Brochard and Nguyen's work [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF]. Furthermore, the kinematic, dynamic and thermal relations at the arc-pool interface are carefully derived and the dependence of the dynamic and thermal equations onto the kinematic condition is briefly discussed. Moreover, both the fixed and deformable arc-pool interface hypotheses are presented.

In chapter 3, the numerical methods, coupling techniques and the algorithms used to solve the physical model are presented. A brief discussion of the numerical tools, multiphysics methodologies and the cou-pling algorithms generally available for the resolution of non-linear models is presented. Next, a detailed description of a monolithic mixed heat equation approach is discussed and the stability and conservative nature of the method is shown. Subsequently, a detailed description of the Dirichlet-Neumann and the quasi-monolithic mass and momentum transfer coupling algorithms is discussed. Similarly, the stability and conservatiive nature of the two methods are briefly discussed.

Next, in chapters 4 and 5, 2D axisymmetric configurations are set up for the analysis of the TIG welding model. At affordable computational costs, the sensitivity of the model to both physical parameters and coupling algorithms are studied. Chapter 4 is based on a parametric study that analyses the strong influence that the cathode tip geometry has onto the arc-plasma. To better interpret the dominant physics at play, the anode is maintained solid, while the cathode geometry is modified and its effect on the Lorentz force of the arc is studied. In fact, it is shown that the tip geometry alone can significantly influence the transferred heat flux and imparted arc shear and pressure to the anode. Between the bluntest and sharpest cathodes studied, a 1.7-fold and 4-fold rise in the heat flux and arc stress was calculated, respectively. Thus, this study highlights the importance of incorporating the cathode into the calculation domain, which is not common in much of the literature on TIG modelling. Chapter 5 is based on multiple studies that focus on both the numerical and physical consequences of the arc-pool coupling problem. Hence, an algorithm comparison study is performed where the Dirichlet-Neumann and the quasi-monolithic algorithms are shown to converge to the same solution (to a small order of precision), and are thus considered comparable. For the cases considered, the Dirichlet-Neumann algorithm is shown to outperform the quasi-monolithic algorithm by at least 20% and at most 40%, depending on the relaxation factor used. Thus, the Dirichlet-Neumann algorithm was retained for use in the rest of the work. Furthermore, the influence of the inlet electric current, the interface hypothesis and the chosen liquid steel viscosity is studied. Effectively, the arc-plasma is found to be greatly influenced by inlet electric current, which then significantly affects the geometry and dimensions of the weld pool. Moreover, even though the Weber number of the studied configurations is generally low, the interface hypothesis is found to influence the depth of the weld pool. Lastly, the choice of liquid steel viscosity in the model is shown to significantly affect the shape and depth of the weld pool.

In fact, the use of an augmented pool viscosity value, as considered by many authors [START_REF] Kim | Heat and fluid flow in pulsed current gta weld pool[END_REF][START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF][START_REF] Traidia | Multiphysics modelling and numerical simulation of GTA weld pools[END_REF][START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF], was shown to dampen secondary circulation zones which strongly affect the penetration and to a lesser extent the widening of the pool.

Lastly, two 3D calculations are presented in chapter 6. The chapter begins with a verification study that ensures the coherent extension of the 2D axisymmetric model to a 3D cartesian one. Effectively, a 3D spot configuration is shown to maintain the rotational symmetry of the spot weld while also agreeing with a reference 2D axisymmetric calculation. This approach allowed to verify all implemented and modelled phenomena, including but not limited to, the deformation of the fluid interface. Using the verified 3D model, a welding configuration with modelled displacement effects is constructed. It is shown that the implemented coupling algorithms remain conservative for such welding configurations and that asymmetry in the arc fields reinforces the importance of a fully coupled 3D arc-pool model in TIG welding. Furthermore, a comparison of the calculated weld pool dimensions to an experimental case studied in Koudadje's work [START_REF] Koudadje | Etude expérimentale et modélisation numérique du bain de fusion en soudage TIG d'aciers[END_REF] is performed. The general agreement of the calculated and measured dimensions of the weld pool indicates the general validity of the developed model.

In conclusion, the work presented in this thesis served to develop, implement, verify, and analyse a fully coupled cathode-arc-pool welding model applicable in both 2D axisymmetric and 3D cartesian configurations. Furthermore, this thesis highlights the importance of a fully coupled 3D model for the simulation of TIG welding with torch displacement effects.

Perspectives

The experience gained while working on this thesis helps identify certain perspectives that could further enhance and extend the model developed in this work. The suggested perspectives are thus listed as follows:

1. An extensive validation study, set up in an analogue manner to the 3D weld displacement model, is strongly suggested. Controlling the displacement of the workpiece w.r.t the torch, and as many other welding parameters as possible, should allow to study the extent of the validity of the model.

Moreover, the validation study could offer more insight on any potential limitations of the model.

2. The elevated Reynolds numbers calculated in weld pools encourage deeper investigations of their hydrodynamic regimes. In fact, the work of Kidess et al. [START_REF] Kidess | The influence of surfactants on thermocapillary flow instabilities in low prandtl melting pools[END_REF] shows that certain weld pools can develop turbulent characteristics. To capture the various thermohydraulic phenomena, they implement a Large Eddy Simulation model in the pool domain. Thus, the investigation and consideration of some appropriate turbulent model for the weld pool is encouraged.

A.2 Temperature dependent material properties

Here the relevant thermophysical of the four used materials in this thesis are presented as they are assembled and/or reported by different authors. The temperature dependent electric and thermal properties of argon plasma, thoriated tungsten (2% Th), solid copper and 316L steel are plotted and presented. Moreover, the viscosity of liquid steel is presented as reported and/or used by different authors.

Argon plasma

The electrical and thermal conductivities, the mass density, viscosity, specific heat capacity and net plasma emission coefficient of argon in its gas and plasma states as functions of temperature. •sr -1 ], respectively. Data is assembled in Brochard's thesis [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF].

Thoriated tungsten

The electrical and thermal conductivities of solid thoriated tungsten as functions of temperature. ] conductivities of each respective material. Data is assembled in Brochard's thesis [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF].

Solid copper

The electrical and thermal conductivities of solid solid copper as functions of temperature. , where the temperature is the primal or dependent variable [START_REF] Nedjar | An enthalpy-based finite element method for nonlinear heat problems involving phase change[END_REF][START_REF] Vaughan R Voller | An enthalpy method for convection/diffusion phase change[END_REF]. A more detailed discussion of isothermal phase change is out of the scope of this text, and so the interested reader is referred to Voller [START_REF] Vaughan R Voller | An enthalpy method for convection/diffusion phase change[END_REF], Nedjar [START_REF] Nedjar | An enthalpy-based finite element method for nonlinear heat problems involving phase change[END_REF], Kumar [START_REF] Kumar | 5.12 -casting simulation methods[END_REF] and Ayasoufi [3].

A discussion of isothermal phase change sets the stage for a discussion of phase change as it occurs for alloys. This thesis being mainly concerned with steel alloys, the mushy zone at the phase transition zones nature of alloys (see [START_REF] Kurz | Fundamentals of Solidification: Fourth Revised Edition[END_REF][START_REF] Kumar | 5.12 -casting simulation methods[END_REF]). Thus, analogously to equation (B.2), the enthalpy-temperature relation can be modelled by the following piecewise function: It is apparent that the presence of a mushy zone in the phase change process should render the problem less stiff when solving it numerically [START_REF] Nedjar | An enthalpy-based finite element method for nonlinear heat problems involving phase change[END_REF]. In fact, unlike isothermal phase change, the phase transition mushy zone manifests as a finite function representing the transition phase 3 ∂L f ∂T in the enthalpy-temperature relation. The transition phase is consequently modelled as: and from the theory of distributions (see [START_REF] Robert | A guide to distribution theory and Fourier transforms[END_REF][START_REF] Piela | Appendix E: Dirac Delta Function -Ideas of quantum chemistry[END_REF] for details), a normalised Gaussian distribution behaves as a δ-Dirac in the limit that its width σ → 0, as described below:

h(T ) =                      T T ref c
∀ L f = ∞ -∞ ∂L f ∂T dT ∃ ∂L f ∂T = L f f (T ; ∆T m , T c ) then ∞ -∞ ∂L f ∂T dT = L f
δ(T ) = lim

6 CHAPTER 1 .

 61 Figure 1.1: Principal constituents of the TIG welding system.

I

  begin by introducing the scaling arguments concerned with the tangential force balance at the arc-pool interface: Tangential interface dynamics: (τ pool ) * = SAP (τ arc ) * + M P S(∇ s T ) * (1.1) where the * superscript represents non-dimensionality, τ the tangential viscous stresses at the arc-pool interface, SAP the shear arc-pool and M P S the Marangoni pool shear numbers, respectively. The normal force balance at the interface are: Normal interface dynamics: (∇ • n) * = W e pool (P pool ) * -W e arc (P arc ) * -Bo(∆ρ) * + Ca pool (σ pool ) * -Ca arc (σ arc ) * (1.2)

Figure 1 . 2 :

 12 Figure 1.2: TIG spot configuration representing the dynamic fluid interface.

Figure 2 . 1 :

 21 Figure 2.1: Schematic of the geometry of the TIG system.

Figure 2 . 2 :

 22 Figure 2.2: Gaussian pillbox over the arc and pool.

Figure 2 . 3 :

 23 Figure 2.3: Level set at the arc-pool interface.

CHAPTER 2 .Figure 2 . 4 :Figure 2 . 5 :

 22425 Figure 2.4: Sum of forces at a fixed interface element.

Figure 2 . 6 :

 26 Figure 2.6: Superimposed plots of schematised pool specific heat c p and enthalpy h. Not drawn to scale.

Contents 3 . 1 42 3. 2 42 3. 3 A 46 3. 4 49 3. 4 . 2 51 3. 5 53 3. 5 . 2 57 3. 6 59 3. 7 66 3. 8 68 3. 1 Introduction

 3142242346449425155352576597668681 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Geometric preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . brief discussion of the chosen numerical methods . . . . . . . . . . . . . . . . 42 3.3.1 The Galerkin FEM approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3.2 The LSFEM approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The multiphysics computational toolbox . . . . . . . . . . . . . . . . . . . . . . 48 3.4.1 Coupling techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Solution algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conjugate heat transfer coupling algorithm . . . . . . . . . . . . . . . . . . . . 53 3.5.1 Interfacial thermal conditions, a mixed heat equation approach . . . . . . . . . . . Identifying the pool boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.5.3 A note on the stability of the interfacial temperature . . . . . . . . . . . . . . . . . Momentum and mass transfer coupling algorithm . . . . . . . . . . . . . . . . . 58 3.6.1 The algebraic systems of equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.6.2 Coupling schemes at the fluid interface . . . . . . . . . . . . . . . . . . . . . . . . . Convergence criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . This chapter discusses the numerical methods, coupling techniques and algorithms as they are implemented and constructed for use in this work. The chapter begins with geometric preliminaries, discussing the numerical construction of certain aspects in the TIG system; and follows with section 3.3 which briefly presents the different finite element methods (FEM) used to treat the different modelled equations in this work; then follows with section 3.4 which briefly presents the coupling techniques and solution algorithms used to build the multiphysics algorithm implemented in this work. Next, a discussion of the construction of a mixed heat equation algorithm is presented in section 3.5.1, where the temperature based heat equation and the enthalpy based one are coupled. Section 3.6 discusses the implementation of a Dirichlet-Neumann algorithm to solve the rigid, partitioned fluid system, and a quasi-monolithic algorithm used to solve the rigid fluid system. Moreoever, the treatment of the deformation at the fluid interface is discussed. Before concluding the chapter, the criteria used to asses the convergence of a stationary solution are presented, as they are defined for both 2D and 3D simulations. The chapter is concluded with a flowchart, representing the multiphysics algorithm, coupling all physics modules, as used in this work.

Figure 3 . 1 :

 31 Figure 3.1: Simplified schematic of constituents of global domain Ω T IG .

4 .

 4 As for the boundary counditions of the problem: a) the Dirichlet type conditions are applied as equivalent fluxes in the {F } term at the nodes across the ∂Ω D of the mesh using the Lagrange Multiplier technique 1 [64], and b) Neumann type conditions are applied as surfacic fluxes in the {F } term at the nodes pertaining to the ∂Ω N of the mesh. First, N ψ ∈ H 1 0 , which means that the interpolation functions N ψ are square integrable and are null at boundaries N ψ | Γ = 0. The form of equation (3.1) as it is set up for the electric, heat transfer and fluid dynamics models is thus presented in the following:1. The electric model:

  {F } are the Neumann boundary conditions imposed as equivalent fluxes. System (3.2) is weakly non-linear 2 , and so it requires coupling to the thermal equation and linearisation before solving for {ϕ}. Note that in the Ω pla , the σ * in the direct vicinity of Γ CP I and Γ AP I is assigned the conductivity values as they are calculated at a δl cat = 0.1 mm and δl ano = 0.4 mm Eucledian distance from the respective interfaces, regardless of the temperature field calculated at T (r < δl cat,ano ) andof the meshing in that zone. This method is adapted from Lowke[START_REF] Lowke | LTE diffusion approximation for arc calculations[END_REF] and is schematised in figure3.2, representing how σ * is calculated in the vicinity of the Γ CP I . The same mathematical operation is applied in the plasma in the vicinity of the Γ AP I interface.

Figure 3 . 2 :

 32 Figure 3.2: Visualisation of how σ * is calculated in the plasma in the vicinity of the Γ CP I interface.

Figure 3 . 3 :

 33 Figure 3.3: Schematic of constructed modules of the global model. Each module transmits physically interpretable data in the form of velocities, temperatures, forces, fluxes etc....

Figure 3 .

 3 4 schematises the partitioning by domain of the fluid dynamics module. The communication of velocities and forces between the submodules occurs at the partition interface, which for this module, occurs at the Γ AP I interface.

Figure 3 . 4 :

 34 Figure 3.4: Fluid dynamics module further partitioned by domain. The arc and pool submodules are connected by interfacial conditions, and after resolution are connected to the interface deformation submodule.

  [START_REF] Carlos A Felippa | Partitioned analysis of coupled mechanical systems[END_REF] where [T hh ] k-1 represents the enthalpy tangent matrix, {δh} k and {R h } k-1 the enthalpy and residual vectors, defined across the nodes of the meshed workpiece subdomain, respectively. The tangent matrices of linear systems (3.24) and (3.25) are detailed in appendix C. Linear system (3.24) contains the energy conservation equations for both the cathode and the arc plasma. Linear system (3.25) expresses the energy conservation equation as it manifests in its enthalpic form and treats only the workpiece domain. The boundary conditions of the enthalpy model are assumed included in linear system(3.25). As for the arc-pool interface Γ AP I , it is treated using the Lagrange Multiplier technique implemented in the Cast εM toolbox[START_REF]Cast3m Web site[END_REF][START_REF] Pellet | Dualisation des conditions aux limites[END_REF]. Thus, the communication between system (3.24) and(3.25) is what is left to treat before solving for the primal fields {T } k and {h} k . The communication occurs at the interface between the arc and the workpiece.

Figure 3 .Figure 3 . 5 :

 335 Figure 3.5: Section of Γ AP I interface where interfacial conditions are applied. Domain i represents arc nodes at Γ plaapi , and j anode nodes at Γ anoapi and are made up of N interfacial nodes.

  [START_REF] Gonzalez | Numerical modelling of an electric arc and its interaction with the anode: part II. the three-dimensional model-influence of external forces on the arc column[END_REF] where system (3.29) is a well-posed linear system of equations. Note that the interfacial equations(3.27) and(3.28) are only applied across spatially overlapping nodes. The use of condition (3.27) and constraint (3.28) strongly 9 couples linear systems (3.24) and (3.25) [88]; however, regardless of the monolithic scheme, the use of different primal variables does not rigidly ensure the continuity of temperature at the interface as required by equation (2.74). Rather, the coupling relations ensure the continuity of equation (3.27) at iteration k. This loosely imposed temperature continuity is considered satisfactory because as the global algorithm converges, relation (3.27) behaves like:

Figure 3 . 6 :

 36 Figure 3.6: Typical temperature-enthalpy relation of an non-isothermal phase change process.

Figure 3 . 7 :

 37 Figure 3.7: (a) Maximum temperature at Γ plaiap (black plot) and at Γ anoiap (red plot) interfaces as functions of iteration number. (b) Convergence plot of all variables of sample simulation.

Figure 3 . 8 :

 38 Figure 3.8: Section of Γ AP I interface where Dirichlet-Neumann coupling is applied. Domain i represents arc nodes at Γ plaapi , and j anode nodes at Γ anoapi and are made up of N interfacial nodes.

35 )Figure 3 . 9 :

 3539 Figure 3.9: Domain i represents arc nodes, and j anode nodes for all N nodes at the interface. Blue arrows represent translation of calculated α k i forces from Γ plaapi to Γ anoapi . Red arrows represent inverse translation for velocities v k j .

3 . 40 ) 3 .Figure 3 . 10 :

 3403310 Figure 3.10: Steps required to solve the Dirichlet-Neumann coupling algorithm. Superscript m is the subcycle iteration counter.

Figure 3 . 11 :

 311 Figure 3.11: Maximum velocities at Γ plaiap (black plot) and at Γ anoiap (red plot) interfaces as functions of iteration number.

Figure 3 . 12 :

 312 Figure 3.12: Red mesh Ω cat , green Ω pla , blue Ω ano and red contour ∂Ω pool . (a) Rigid mesh with the melt pool contour. (b) Deformed mesh with melt pool contour.

Figure 3 . 13 :

 313 Figure 3.13: Algorithm of global model. Each block represents a physics module, and the modules developed in this thesis are schematised in further detail. Dashed arrows represent data flow from iteration k -1, and the solid arrows, data from iteration k.

  between the segments [AZ], [ZG], [GF], with tip truncation angle α, radius at the tip r int . The range of the parameters are: α = [15 o , 20 o , 25 o , 30 o ] and r int = [0.15 mm, 0.3 mm, 0.6 mm].

Figure 4 .

 4 1, b) presents a superimposition of the three considered cathode shapes. The axis of symmetry in figure 4.1 a) and b) is indicated by the dashed lines with dots. The parameterisation of the geometry yields a total of 28 cases for the present parametric study. The widths and heights of the computational subdomains Ω cat , Ω pla and Ω ano are kept constant throughout the study. The geometric parameters of the cathode are presented as follows (cf. figure 4.1):

table 4 . 1 ,Figure 4 . 1 :

 4141 Figure 4.1: Sketches of the geometric parametric study. a) The computational domain. b) A close-up view of the three considered cathode shapes, superimposed.

j

  imp • n dΓ EF = I imp with I imp = 200 A → j imp = I imp Γ EF • n and where the electric problem is approached as a fixed current circuit (at 200 A), from which the voltage is calculated. Regarding the thermal boundary conditions, the zero diffusive flux condition is used to apply both a symmetry condition along boundary [EH] and an exit condition along boundary [DC]. The imposed temperature is chosen based on that used by different authors [33, 30].

  σ * , λ, etc...) → Q1, P1 (ϕ, T, B, u) → Q2, P2

Figure 4 . 2 :

 42 Figure 4.2: The meshes for the PNT α = 15 o case. a) M 1 : coarsest mesh with 15, 487 elements b) M 3 : finest mesh with 42, 089 elements. The red, green and blue zones represent the meshes of the cathode, arc-plasma and anode, respectively.

4 )Figure 4 . 3 :

 443 Figure 4.3: Velocity profiles for the three meshes considered and various cathode shapes. Legend: the first number indicates mesh number, the second one the truncation angle (in degrees) and the third one the tip radius [mm]. The axis origin is at the API.

  , so any modification to this force translates into the behaviour of the arc. The geometry of the cathode tips strongly influences the current density at the CPI, as shown in figure4.4, which presents the electric current density j CP I along the curvilinear abscissa at the CPI. The sharper the cathode tip, the more intense the current density profile near the tip, with the maximum current density value rising more than one order of magnitude as compared to the other cases. The different current densities along the CPI result in strongly different

Figure 4 . 4 :

 44 Figure 4.4: The calculated j CP I for the 3 sharpest cathode cases vs. the curvilinear abscissa along the CPI, where the origin is placed at the tip. The legend lists the geometry type, truncation angle and tip radius, respectively.

Figure 4 . 5 :Figure 4 . 6 :

 4546 Figure 4.5: Contour plots of norms of Lorentz force (top row) and velocity field (bottom row) around two different cathode tips. PNT α = 15 o (left column) and RND α = 20 o , r int = 0.6 mm (right column) are presented. Lorentz force scale is in [N•m -3 ] and velocity field in [m•s -1 ].

Figure 4 . 7 :

 47 Figure 4.7: Heat transfer modes at the API for the RND α = 30 o , r int = 0.6 mm case. Negative values represent the heat flux leaving the anode, and the positive one, entering it.

Figure 4 . 8 :

 48 Figure 4.8: Heat flux transferred to workpiece across the API versus maximum Peclet number in the arc.The color key groups cathodes of same angle, while the marker key groups cathode shape and tip size.

Figure 4 . 9 :

 49 Figure 4.9: Maximum arc velocity versus truncation angle (in degrees) and tip type.

  8-fold respectively. The total heat flux into the anode increases ≈ 1.2-fold. These trends are visualized in figure 4.10. The numerical profiles and the fitted functions to the normalized pressure, shear stress and heat flux variables are presented for the PNT α = 15 o and the CHF α = 30 o r int = 0.6 mm cases in figures 4.10 a), b) and c), respectively. The integrated pressure and shear forces are plotted in figures 4.10 d), e) and the integrated heat flux into the anode in f). The maxima and the fitted widths of the pressure, shear and heat flux values of all the simulated cases in this study are presented in figure 4.11. Figures 4.11 a) and b) present the maxima and Gaussian variance of the fitted pressure profiles; figures 4.11 c) and d) the maxima and Weibullian scale of the fitted shear profiles; figures 4.11 e) and f) the maxima and Gaussian variance of the fitted heat flux profiles respectively, of all simulated cases in this study. The shape parameter k is approximately constant for all Weibull distributions fitted to the shear profiles, with k ≈ 1.4. The variability of the electric current density profile j AP I at the API is calculated to be weak for a 10 mm arc. The j AP I profiles across different cases is negligible where the maximum current densitiy j max and the fitted Gaussian widths σ j are seen below: 2.85 × 10 6 ≤ j max ≤ 3.11 × 10 6 [A • m -2 ] (4.7) 3.20 ≤ σ j ≤ 3.34 [mm] (4.8)

Figure 4 .Figure 4 .

 44 Figure 4.10: a) Normalized Gaussian fits and pressure profiles; b) normalized Weibull and shear stress profiles; c) normalized Gaussian fits and heat flux profiles for PNT α = 15 o and CHF α = 30 o , r int = 0.6 mm cases. Black solid lines and red dashed lines indicate the fits calculated profiles, respectively. Axes are non-dimensionalized w.r.t PNT α = 15 o case. d) The integrated pressure force; e) the integrated viscous shear force; f) the integrated heat flux versus cathode geometry.

4. 11 ,

 11 where their results overlap. However, although the electric current density locally differs between the CHF and RND (cf. figure 4.4) this localized difference does not dominate the entirety of the profile, and only makes up a couple of Amperes in electric current. Therefore, the Lorentz force resulting from the CHF and RND cases are very similar, cf. figure 4.12, which presents the contour plots of the Lorentz force in the vicinity of the RND and CHF cases at r int = 0.3 mm, α = 15 o .

Figure 4 . 12 :

 412 Figure 4.12: Contour plots of the norms of the Lorentz force around two different cathode tips. The RND (right) and CHF (left) α = 15 o , r int = 0.3 mm cases are presented. The scale is in [N•m -3 ].

Figure 4 . 13 :

 413 Figure 4.13: Mechanical power of the Lorentz force in the arc versus cathode geometry. Each black shape groups cases that have quasi-equivalent power values.

Figure 4 .

 4 14 presents the voltage drop across the entire domain versus cathode geometry. The voltage drop increases w.r.t increasing tip sharpness. This proportionally increases the total energy available in the circuit as P circuit = I∆ϕ [W], with ∆ϕ the voltage drop, and I the input current at a constant 200 A. The voltage drop as a function of α and r int is presented in figure 4.14. The more the electric power is pulled by the circuit, the more the heat is deposited at the anode. Correspondingly, figure 4.10 (d) shows that the sharper the cathode tips, the larger Q ano is at the work piece. The ratio of the rise of energy deposited by the arc and power pulled into the circuit is quantified by calculating the efficiency η of the circuit. The equation used to calculate the circuit efficiency takes into account the source and sink terms injected at the interfaces at both the CPI and API is:

Figure 4 .

 4 15 a) sketches the electric current density distribution as imposed at the cathode tip in Goodarzi's work, and the one resulting from the computation from the present model.The magnitudes of the main variables of the TIG arc system are presented in table 4.4 for both the Goodarzi et al. study[START_REF] Goodarzi | The effect of the cathode tip angle on the GTAW arc and weld pool: I. mathematical model of the arc[END_REF] and ours. Comparing the maximum heat flux pulled into the anode between the bluntest and sharpest cases, Goodarzi et al. calculated negligible variation (≈ 0.15%), whereas the computed results reveal a ≈ 1.7-fold increase. Moreover, referring to the maximum arc velocities, pressures and shear stresses, the values calculated by Goodarzi et al. vary significantly less than ours. Furthermore, I note that the electric current density, at the cathode tip, is quasi-constant for the range of cases studied by Goodarzi et al.. Indeed, as previously stated, the electric current density in the vicinity of the cathode tip directly influences the Lorentz force. Therefore, a lack of variability in the electric current density1 They assign values to the cone angle, whereas in the present work, the half angle is used.

Figure 4 .

 4 Figure 4.15: a) Schematic representations of electric current emission zone for Goodarzi's cathode, and ours. b) Estimated spot areas in this study, those used in Goodarzi [30] and those from Haidar [32] versus cathode geometry. The error bars from Haidar's work are not displayed.

  behave in a purely diffusive manner. The studied configurations consider a 10 mm arc height, a 200 A inlet electric current, for three tip shapes: pointed, chamfered and rounded tips. The range of truncation angles and radii are: α = [15 o , 20 o , 25 o , 30 o ]; and r int = [0.15, 0.3, 0.6] mm. The main findings of this section are three-fold. First, the shape of the cathode tip strongly influences the current density at the cathodeplasma interface, that in turn significantly affects the resulting Lorentz force. This being the main source of momentum in the fluid flow of the arc-plasma, it induces a broad variability of the transported quantities relevant to the welding parameters. Indeed the heat flux pulled in by the anode shows a direct dependence on the Peclet number, which reflects the magnitude of the fluid flow in the arc-plasma. The heat flux varies up to 1.7 times, while the pressure and shear stresses rise about 4 times, respectively, from the bluntest to the sharpest tips. Therefore, this parametric study strongly suggests the importance of considering

  symmetric configuration. In referring to the figure the rotational axis of symmetry is represented by the dashed line. The cathode truncation angle and radius, and the arc height are maintained constant at α = 15 o , r int = 0.3 mm and h arc = 5 mm, respectively. The nozzle wall is emulated by including the segment [JK] into the domain, where zero flow is allowed through. The inlet flow is imposed between the nozzle wall and the cathode domain at segment [FJ].

Figure 5 . 1 :

 51 Figure 5.1: Schematic of the geometry used.

2 )

 2 Viscosity extracted from the IAEA report[START_REF] Bobkov | Thermophysical properties of materials for nuclear engineering: a tutorial and collection of data[END_REF]

Figure 5 . 2 :

 52 Figure 5.2: Global mesh with a superposed schematic of expected pool zone at successive refinements.

  3 mm c S = 10 ppm µ → equation (5.1) Interface hypothesis → deformable Case 2:

Figure 5 . 3 :

 53 Figure 5.3: Comparing meshes for case 1. Top figure is the pool geometries of the three meshes superposed and is drawn to scale. Bottom figures are, from left to right, temperature at API, velocity norm at API and vertical velocity component at symmetry axis, respectively.

Figure 5 . 4 :

 54 Figure 5.4: Comparing meshes for case 2. Top figure is pool geometries of the three meshes superposed and is drawn to scale. Bottom figures are, from left to right, temperature at the API, velocity norm at API and vertical velocity component at symmetry axis, respectively.

3 .

 3 The fourth setup is based on the I = 125 A configuration run using a relaxation factor of ω N S = 1.0 on both algorithms. Summarising, the total number of cases are named and presented in the table below, where D-N and Q-M stand for the Dirichlet-Neumann and quasi-monolithic algorithms, respectively: Before running a simulation, each case is loaded with an initial condition 1 where the pool hydrodynamics are turned off. For example, cases 1.a and 1.b are initialised with the same solution to a calculation at I = 100 A and for which the anode is maintained ω N S = 0.3 ω N S = 1.0 I = 100 I = 125 I = 150 I = 125 D-N Case 1.a Case 2.a Case 3.a Case 4.a Q-M Case 1.b Case 2.b Case 3.b Case 4.b Table 5.3: Summary of the eight cases used in comparing algorithm performance.

Case 1 :Figure 5 . 5 :

 155 Figure 5.5: a) Convergence plot of case 1.a); b) convergence plot of case 1.b). Legends indicate the increments of primal variables per iteration number. c) Profile of velocity norm at API; d) profile of zcomponent of velocity at the pool symmetry axis.

Case 2 :Figure 5 . 6 :

 256 Figure 5.6: a) Convergence plot of case 2.a); b) convergence plot of case 2.b). Legends indicate the increments of primal variables per iteration number. c) Profile of velocity norm at API; d) profile of zcomponent of velocity at the pool symmetry axis.

Case 3 :Figure 5 . 7 :

 357 Figure 5.7: a) Convergence plot of case 3.a); b) convergence plot of case 3.b). Legends indicate the increments of primal variables per iteration number. c) Profile of velocity norm at API; d) profile of zcomponent of velocity at the pool symmetry axis.

Case 4 :Figure 5

 45 Figure 5.8: a) Convergence plot of case 4.a); b) convergence plot of case 4.b). Legends indicate the increments of primal variables per iteration number.

Case 4 .

 4 b) converges for a smaller number of global configurations when compared to the analogous cases 2.a) and 2.b). However, the total time to convergence of case 4.b) is longer than those of cases 2.a) and 2.b). The solution at convergence for case 4.b) is consistent with the solutions found for the analogous cases 2.a) and 2.b).

  [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF] [START_REF] Sheldon | Methodology for comparing coupling algorithms for fluid-structure interaction problems[END_REF], the relaxed Dirichlet-Neumann algorithm is the most adapted to the work in this thesis. While the quasi-monolithic algorithm in this work can require a smaller number of iterations until convergence w.r.t the Dirichlet-Neumann algorithm, it performs slower.

5. 6 . 2 . 1 Figure 5 . 9 :Figure 5 .

 621595 Figure 5.9: a) z-component of velocity along the arc symmetry axis. b) Temperature profile along the arc symmetry axis.

figure 5 .

 5 figure 5.13 and the maps of the stream function field of the I = 100 A and I = 150 A cases (figures 5.12.a) and b)), the following arguments about the influence of electric current onto pool width and depth are made.

Figure 5 .

 5 Figure 5.11: a) Velocity norm profiles at API. b) Z-component of velocity at the pool symmetry axis. c) Temperature profiles at API. d) Temperature profiles at the pool symmetry axis.

Figure 5 .

 5 Figure 5.12: a) Stream function field with superimposed velocity vectors for the I = 100 A simulation. b) Stream function field with superposed velocity vectors for the I = 150 A simulation.

Figure 5 . 13 :

 513 Figure 5.13: Pool geometries of the three simulations run for different inlet electric currents.

Figure 5 . 14 :

 514 Figure 5.14: Marangoni surfacic force at pool surface.

Case 1 :

 1 I arc = 150 A α = 15 o r int = 0.3 mm c S = 10 ppm µ → equation (5.2) Interface hypothesis → fixed Case 2: I arc = 150 A α = 15 o r int = 0.3 mm c S = 10 ppm µ → equation (5.2) Interface hypothesis → deformable.

Figure 5 .

 5 Figure 5.15: a) z-component of velocity at the pool symmetry axis. b) Pool geometries of the I = 150 A cases with and without the Lorentz force.

r

  int = 0.3 mm c S = 10 ppm µ → equation (5.1) (augmented viscosity) Interface hypothesis → deformable Case 2: I arc = 150 A α = 15 o r int = 0.3 mm c S = 10 ppm µ → equation (5.2) (experimental viscosity) Interface hypothesis → deformable.

Figure 5 . 16 :

 516 Figure 5.16: Pool geometries of the two simulations run for different interface hypotheses.

Figure 5 . 17 :

 517 Figure 5.17: Streamfunctions in the region of the secondary circulation zones of case 1 (left) and case 2 (right).

Figure 5 .

 5 Figure 5.18: a) Z-component of velocity along the arc symmetry axis. b) Temperature profile along the arc symmetry axis.

Figure 5 . 19 :

 519 Figure 5.19: Pool geometries of the two simulations run for different liquid steel viscosities.

Figure 5 .

 5 Figure 5.20: a) Velocity norm profiles at API; b) z-component of velocity at the pool symmetry axis; c) temperature profiles at API; d) temperature profiles at the pool symmetry axis.

Figure 5 . 21 :

 521 Figure 5.21: Reynolds number profiles at the pool surface for case 1 and 2.

figure 6 . 1 .

 61 figure6.1. The interface between the arc-plasma and the weld pool, defined along segment [QN], is treated using the Dirichlet-Neumann algorithm (see section 3.6.2.1 in chapter 3).

6. 1 . 4 . 1 Figure 6 . 3 :

 14163 Figure 6.3: Temperature field of the side view of all TIG subdomains (left), side view of weld pool (top right) and top view of weld pool (bottom right).

Figure 6 . 4 :Figure 6 . 5 :

 6465 Figure 6.4: Velocity norm field of the side view of weld pool (top) and top view of weld pool (right).

Figure 6 . 6 :

 66 Figure 6.6: Comparing the two meshes in 3D to the analogous mesh in 2D. a) Temperature profile along arc symmetry axis. b) Vertical velocity component profile along arc symmetry axis. c) Temperature profile along API (center to edge). d) Velocity norm profile along API (center to edge). e) Temperature profile along anode symmetry axis. f) Vertical velocity component profile along pool symmetry axis.

Figure 6 . 7 :

 67 Figure 6.7: Pool contours of mesh 1 and 2 of the 3D model and axisymmetric mesh of the 2D model.

Figure 6 . 8 :

 68 Figure 6.8: Half geometry used in section 6.2, sliced at the z -x symmetry plane, where y = 0.

  implications of the thermohydraulic fields onto a) the algorithms implemented in this work and b) the weld pool geometry. First the temperature field of the global domain is presented in figure 6.11.

Figure 6 . 11 :

 611 Figure 6.11: Temperature field of simulated case of the weld displacement configuration.

Figure 6 . 12 :

 612 Figure 6.12: Temperature field of top view of the the workpiece at Γ AP I .

6. 14 ,Figure 6 . 13 :

 14613 Figure 6.13: Arc and pool temperature and velocity fields. U represents the streampath scale, and V the vector field scale.

Figure 6 . 14 :

 614 Figure 6.14: Pool geometry and the temperature field in its vicinity.

Figure 6 .

 6 Figure 6.15: A closeup of the pool temperature and velocity fields.

Figure 6 . 16 :

 616 Figure 6.16: Calculated maximum depth isocontour on the simulated pool geometry.

3 .

 3 The consideration of the effects of filler metal during TIG welding should be investigated so as to further extend the applicability of the current model towards more industrial configurations. A careful consideration of both the physical and numerical effects of the momentum and mass source terms at the fluid interfaces is necessary. 4. The consideration of transient phenomena in the physical model is encouraged. The transient effects should be developed into the numerical model by the use of implicit integration schemes for their stability. Furthermore, I encourage research on automated relaxation methods (see Aitken's deltasquared method as an example) which could potentially accelerate and/or enhance the stability of the numerical model [49, 23, 20] . 5. Finally, the development of remeshing strategies for use in the vicinity of the arc-pool interface and in the weld pool is encouraged. Remeshing would potentially allow for quicker numerical convergence of the fully coupled model, especially when configurations with significant fluid interface deformation are simulated.

Figure A. 1 :

 1 Figure A.1: Material properties as functions of temperature T [K]. Subfigures (a) and (b) present the electrical σ [S•m -1 ] and thermal λ [W•m -1 K -1 ] conductivities, respectively. Subfigures (c), (d), (e) and (f) present the density ρ[kg•m -3 ], dynamic viscosity µ [kg•m -1 •s -1 ], heat capacity c p [J•kg -1 •K -1 ] and the net plasma emission coefficient ϵ n [W•m -3 •sr -1 ], respectively. Data is assembled in Brochard's thesis[START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF].

Figure A. 2 :

 2 Figure A.2: Material properties as functions of temperature T [K]. Subfigures (a) and (b) tungsten present the electrical σ [S•m -1 ] and thermal λ[W•m -1 •K -1 ]conductivities of each respective material. Data is assembled in Brochard's thesis[START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF].

Figure A. 3 :

 3 Figure A.3: Material properties as functions of temperature T [K]. Subfigures (a) and (b) present the electrical σ [S•m -1 ] and thermal λ[W•m -1 •K -1 ]conductivities, respectively. Data is assembled in Brochard's thesis[START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF].

Figure B. 1 :T(B. 4 )

 14 Figure B.1: Typical enthalpy-temperature relation of an isothermal phase change process.

Figure B. 2 :

 2 Figure B.2: Typical specific heat-temperature relation of an isothermal phase change process.

  ps dT for T ≤ T s Ts T ref c ps dT + T -∞ ∂L f ∂T dT for T s < T < T l Ts T ref c ps dT + T T l c p l dT + L f for T > T l (B.5)

  ; ∆T m , T c ) dT for ∆T m = T l -T s for T c = T l + T s 2 where in the limit that ∆T m → 0, the enthalpy-temperature relation approaches the isothermal transition relation of a pure metal (see equation (B.2)). Thus: ∆T m , T c ) dT = L f ∞ -∞ δ(T -T c ) dT (B.6)

T c ) 2 2σ 2 )

 2 Figure B.4: Typical specific heat-temperature relation of an non-isothermal phase change process.

process 4 .

 4 Figure B.6: Typical h -T relation schematicising integration over the mushy zone using four mesh nodes.

Table 1 .

 1

	Dimensionless number Symbol Force Ratio	Equation
	Shear arc-pool	SAP	Arc to pool shear	(U a µ a /L a ) (U p µ p /L p )
	Marangoni pool shear M P S	Marangoni to pool shear	| dγ dT |	∆T (U p µ p )
	Arc Weber	W e arc	Arc dynamic pressure to surface tension	P a (γ/L p )
	Pool Weber	W e pool	Pool dynamic pressure to surface tension	P p (γ/L p )
	Bond	Bo	Gravity to surface tension	(ρ p -ρ a )gL p γ/H z
	Arc Capillary	Ca arc	Arc viscous stress to surface tension	U a µ a L p L a γ
	Pool Capillary	Ca pool	Pool viscous stress to surface tension	U p µ p γ

1: Description and equations of non-dimensional numbers at arc-pool interface.

Table 1 . 2

 12

: Magnitudes of characteristic non-dimensional numbers at arc-pool interface. Zones a,b and c refer to figure 1.2

Table 3 .

 3 1: Convergence criteria used in this work.

Table 4 .

 4

1: Applied boundary conditions (cf. figure 4.1).

Table 4 . 2

 42

: Table of used elements and the variables associated to them.

  Three meshes have been generated for each case, a coarse, intermediate and fine denoted as M 1 , M 2 and M 3 , respectively (cf. figure 4.2 that presents M 1 and M 2 of the PNT α = 15

o case). The unstructured meshes are built with both triangular and quadrilateral elements, having 7 and 9 nodes each, respectively.

The high element density zones in the meshes are generated in regions where high gradients are expected, namely close to the cathode tip, symmetry axis and the API. Table

4

.3 shows the number of elements

Table 4 .

 4

	PNT	PNT	CHF, all r int	RND, all r int
	15 o 20 M 1 = 15,487 ∼ 4,000	∼ 3,500	∼ 3,500
	M 2 = 22,786 ∼ 7,000	∼ 7,000	∼ 7,000
	M 3 = 42,089 ∼ 13,500	∼ 14,000	∼ 14,000

o ≤ α ≤ 30 o 15 o ≤ α ≤ 30 o 15 o ≤ α ≤ 30 o 3:

The approximate number of elements of the 3 meshes used for all studied cases.

Table 4 . 4

 44

	Case	Bluntest Sharpest Ratio Bluntest Sharpest Ratio Bluntest Sharpest Ratio
	Goodarzi [30]	≈ 200	≈ 200	1	39.30	39.84	1	398.9	413.8	1.04
	This work	40	30000	750	27	45		1.7	170	470	2.8
		Variable		P max [Pa]			τ max [Pa]
		Case	Bluntest Sharpest Ratio Bluntest Sharpest Ratio
		Goodarzi [30]	850	1078	1.3	125	150	1.2
		This work	267	935		3.5	34	139	4

: Characteristic results from Goodarzi et al.

Table 4 .

 4 5: Results from experimental studies in the literature versus this study. comparison presented in table 4.5 it turns out that the present model is capable of reproducing the general experimental trends concerning the main physical quantities of interest. Moreover, it is interesting to note that Tanaka et al.[START_REF] Tanaka | Anode heat transfer in tig welding and its effect on the cross-sectional area of weld penetration[END_REF] observe, for arc configurations similar to ours, that the electric current density at the API does not significantly vary with cathode tip geometry. Thus, it should be considered that the model implemented in this study is relevant to capturing the effects of cathode tip geometries onto TIG arcs.

	This work	267	935	3.5	34	139	4	27	45	1.7
				Variable		∆ϕ [V]
				Case	Bluntest Sharpest Ratio
			Chihoski [17]	17.9	18.6	1.041
			This work a)	11	13	1.18
			This work b)	18	20	1.11

Table 5 .

 5 2: Table of used elements and variables associated to them.

Table 5 .

 5 5: Table summarising calculation times and total number of iterations to convergence. Ratio of recorded performances are calculated.

		Case 1 Case 2 Case 3	Case 4
	D-N a)				Not Converged
	Calculation Time (s)	23264	30772	55924	162055
	Total Iteration	1021	1158	1749	2000
	Time per Iteration (s)	24	27	32	81
	Avg. Internal D-N Iterations 1	1	1	2.5
	Q-M b)
	Calculation Time (s)	22935	37458	62196	43200
	Total Iteration	998	1388	1932	948
	Time per Iteration (s)	23	27	32	46
	Ratio:
	Time D-N to Q-M	1.014	0.822	0.899	--
	Iteration D-N to Q-M	1.023	0.834	0.905	--

Table 5 .

 5 6: Summary of pool widths, depths and depth to width ratios.

		I = 100 I = 125 I = 150
	Width (mm) 5.95	7.51	9.14
	Depth (mm) 1.76	2.26	3.22
	Ratio D:W	0.30	0.30	0.45

  Geometric configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.1.2 Boundary conditions and material properties . . . . . . . . . . . . . . . . . . . . . 118 6.1.3 Discretisation and meshing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.1.4 Discussion and results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 Geometry configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.2.2 Boundary conditions, displacement effects and meshing . . . . . . . . . . . . . . . 126 6.2.3 Discussion and results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128

	Chapter 6
	A three dimensional study of TIG
	welding
	6.2 TIG welding with displacement effects . . . . . . . . . . . . . . . . . . . . . . . . 125
	6.2.1

Une intelligence qui, à un instant donné, connaîtrait toutes les forces dont la nature est animée et la situation respective des êtres qui la composent, si d'ailleurs elle était suffisamment vaste pour soumettre ces données à l'analyse, embrasserait dans la même formule les mouvements des plus grands corps de l'univers et ceux du plus léger atome ; rien ne serait incertain pour elle, et l'avenir, comme le passé, serait présent à ses yeux. Le démon de Laplace Pierre-Simon Laplace

Contents 6.1 A TIG spot verification study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.1.1 6.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 6.3.1 Spot verification study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 6.3.2 Welding with displacement effects . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

Table 6 . 2

 62

: Used elements and associated variables. Interpolation nodes are represented by the "o" symbol.

The scaling arguments although argued for a TIG spot configuration, are expected to apply for non-stationary configurations.

The fluid interface is deformable if the deformable interface hypothesis is chosen. This is further detailed in the following sections.

This refers to heat concentrated at high temperatures which is advantageous for welding because the temperatures rise significantly above the melting point of the workpiece.

Which ultimately shapes the weld itself.

Sufficiently far from the arc core.

The Gaussian pillbox is a type of Gaussian surface that facilitates the visualisation of Gauss's theorem.

An application of this technique is presented in section 3.6.2.1

Ampere's Law is already a system of three equations and three unknowns. Although this is necessary for the determinacy of an algebraic system, it is not sufficient.

Data means mathematical information, in physically interpretable quantites, i.e: through velocities, temperatures, forces, fluxes etc....

Monolithic at the domain level, meaning that the the field across the entire spatial domain is solved simultaneously.

Allowing for a monolithic approach to the resolution of the energy conservation equations.

The ϵ at the limit is a predefined convergence criterion assumed to be sufficiently small that errors in precision become negligible.

Essentially a type of Gauss-Seidel scheme, see[START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF][START_REF] Patricot | Couplages multi-physiques : évaluation des impacts méthodologiques lors de simulations de couplages neutronique/thermique/mécanique[END_REF]

Before the interface deformation treatment.

The use of the word natural here is to describe a more mathematically consistent approach, where the entire fluid domain is treated as one mathematical system rather than as two, as is done in partitioned methods.

Or iteration m if subcycling occurs.

The term "trial method" is that used in the article this technique is adapted from[START_REF] Cuvelier | Some numerical methods for the computation of capillary free boundaries governed by the navier-stokes equations[END_REF].

What is meant by the "vertical axis" is the unit normal of the initial configuration of the interface, before deformation.

An internal element is an element defined within the domain in question and is different to a boundary element.

The term initial condition is a misnomer here because the model is a stationary one. However, due to the lack of a better term it is used.

The number of degrees of freedom in 3D systems grows much quicker than in 2D systems, and thus the scalability of a direct solver becomes an issue.

A surfacic force is implied here.

Conforming to the mathematical standard, a flow rotating clockwise has a negative stream function associated to it.

The differences between the thermophysical properties of

304L and 316L steel are considered negligible.3 This section was chronologically written last and thus little time was left from my PhD to allocate to expensive calculations.

The integral is equivalent to the latent heat up to a certain order of numerical precision

Acknowledgements

The end of this difficult but ephemeral chapter in my life comes with a plethora of emotions that words cannot describe; thus, my words cannot do justice to the thanks I owe to my family, friends and colleagues. I begin by thanking first and foremost my supervisor Stéphane Gounand who gave me the opportunity to work by his side on this difficult but fruitful project. Times were tough for the both of us, but we came through. I also thank my director Marc Medale for the opportunity to work as a doctoral researcher under

Summary

The different physics modules implemented in this multiphysics TIG welding model serve to simulate the dominant phenomena as they are modelled in this thesis. A mix of numerical techniques and algorithms are used to look for stationary solutions to the global model. Moreover, this chapter proposes a novel coupling algorithm to solve for a mixed heat equation in the global model. The fluid dynamics module is also discussed in detail, where a domain partitioned method and a quasi-monolithic method are discussed and coupled to an interface deformation algorithm. This chapter is concluded with a flowchart summarizing the physics modules as they interact with each other in the algorithm, schematised in figure 3.13, as it is developed to solve the global physical model.

Chapter 4

The influence of cathode geometry on TIG arcs

Understanding is the reward of faith. Therefore, seek not to understand that you may believe, but believe that you may understand.

Saint Augustine of Hippo

Table 6.1: Boundary conditions imposed at the geometry contour of figure 6.1.

As for the material properties used in the verification study, they are also based on a solid thoriated tungsten cathode, argon arc and a 316L steel anode. The thermophysical properties as functions of temperature are presented in appendix A. Note that the chosen viscosity value of liquid steel is based on the augmented viscosity value from the numerical work of Nguyen [START_REF] Chien Nguyen | Sensitivity to welding positions and parameters in GTA welding with a 3D multiphysics numerical model[END_REF]. The viscosity is described by equation (5.1).

Discretisation and meshing

Here the chosen discretizations used for the different variables and the different meshes used are presented.

Note that the discretisation choices presented in the following apply to all discussed simulations in sections

∂Ω anoin ∪ ∂Ω anoout Table 6.3: Boundary conditions imposed at the surfaces of schematised geometry in figure 6.8.

Displacement effects

To account for the weld displacement effects, the workpiece heat transfer model is assigned the displacement velocity w.r.t the laboratory frame of reference V displ everywhere except inside the weld pool domain.

Similarly, the weld displacement velocity is imposed as a moving rigid wall condition (see equation (2.31) in sec. 2.4.3.2, ch.2) where v| ∂Ω pool = V displ . The modelling of the weld displacement effects is schematised in figure 6.9 where the z -x plane of the workpiece at y = 0 is presented. Figure 6.9: Schematic of the modelled displacement effects in the workpiece.

Appendix A

Material properties

A.

Table A.1: Some electric and thermal constants for Tungsten, Argon and Copper. From Brochard's and Pichler's work [START_REF] Brochard | Modele couple cathode-plasma-piece en vue de la simulation du procede de soudage a l'arc TIG[END_REF][START_REF] Pichler | Measurements of thermophysical properties of solid and liquid nist srm 316l stainless steel[END_REF] .

316L steel

The electrical and thermal conductivities, the mass density and the specific enthalpy of 316L in its solid and liquid states as functions of temperature. The liquid steel viscosity as numerically augmented and as it is experimentally reported by Nguyen and an IAEA report, respectively [START_REF] Nguyen | Modelisation et simulation multiphysique du bain de fusion en soudage a l'arc TIG[END_REF][START_REF] Bobkov | Thermophysical properties of materials for nuclear engineering: a tutorial and collection of data[END_REF].

Viscosity extracted from the IAEA report [START_REF] Bobkov | Thermophysical properties of materials for nuclear engineering: a tutorial and collection of data[END_REF] Appendix B

Modelling phase change

Phase change problems in metallurgical applications are generally modelled by non-linear and stiff equations. These equations require interesting workarounds to be applied to their troubling mathematical terms, so as to facilitate their resolution. Therefore, this appendix starts off with a brief discussion of the effects phase change has on to the thermodynamic properties of a given material; and then presents the choice of the energy conservation law adapted to solving the phase change problems seen in the welding of steels, as tackled in this thesis.

Idealized materials and pure metals, under isobaric conditions, can be modelled to undergo phase change across isotherms. These materials undergo what is known as isothermal solidification. The phase change process at constant pressures 1 is characterized by the following equation:

where h is defined as the enthalpy, c p as the specific heat capacity 2 , and T ref as some reference temperature.

Isothermal phase change occurs at a transition temperature, which marks the transformation of the state of the material. Figure B.1 represents a typical phase change relation between the enthalpy of the material and its temperature. For complex materials, such as steels, the enthalpy-temperature relations are typically defined experimentally (see [START_REF] Pichler | Measurements of thermophysical properties of solid and liquid nist srm 316l stainless steel[END_REF] for an example). The discontinuity seen in the enthalpy at the transition temperature T s,l is characterized by a jump in the energy requirement of the process. This jump is due to the higher entropic state of the liquid as compared to that of the solid. To achieve the higher entropic state, energy is spent in liberating the molecular bonds of the material from the low entropy solid state [START_REF] Kurz | Fundamentals of Solidification: Fourth Revised Edition[END_REF]. 1 The discussion is limited to isobaric phase change and thermodynamics as the fluidic systems treated in this thesis are hypothesized to be incompressible. 2 which can be a complex function of temperature not on the same scale as c ps and c p l . The stiffness of the problem would require many integration nodes to be placed in the vicinity of the mushy zone. under the curve represents the numerical approximation of the latent heat L f . Indeed, by looking at the grey numerically integrated region, it can be understood that using the temperature variable in equation (B.9) to precisely capture phase change phenomena is difficult and numerically costly. To overcome the difficulty of integrating over the Gaussian hump, a change of variable and use of the enthalpy as the primal and dependent variable alleviates the numerical stiffness of the phase change problem. This is due to the fact that the enthalpy variable masks the stiffness of the problem because h behaves like an integral quantity, allowing for a smoother problem to be solved. The smoothing of the transition zone is schematised in figure B.6 where four integration nodes are seen to precisely estimate the latent heat contribution. Practically, the change of variable (from T to h) allows for the capturing of the latent heat contribution with modestly sized mesh densities in the phase change zones of a weld pool. Thus, with only a few number of mesh elements in the mushy zone, the latent heat contributions can be precisely estimated.

The change of variable translates equation (B.9) to the following:

As a final note, the peaking behavior of the specific heat is difficult to incorporate into typical numerical heat transfer schemes because the c p term is attached to the transport terms of the conservation equation.

This can introduce large variations in the convection matrix of the system at the mesh nodes in the phase change zone.

Appendix C

Weak formulation

Here the weak formulation of the physical model discussed in chapter 2 is presented. The weak formulation makes up the terms of the numerical model (see chapter 3) that is used to approximate the solutions to the welding configurations studied in this thesis. Moreover, N ψ ∈ H 1 0 , which means that the interpolation functions N ψ are square integrable and are null at boundaries N ψ | Γ = 0.

Nota bene:

If no subscript is indicated in the Ω symbol, the integral is defined across the entire domain. Otherwise, Ω c , Ω p , Ω a represent the cathode, arc-plasma and anode subdomains, respectively.

C.1 Electric model

Non-linear residual

Tangent matrix

C.2 Magnetic model

Residual

Tangent matrix

C.3 Energy conservation model

Non-linear residual

Residual of the temperature based equation, as used in the cathode and arc-plasma subdomains:

Residual of the enthalpy based equation, as used in the anode/workpiece subdomain:

Tangent matrix

Tangent matrix of the temperature based equation:

Tangent matrix of the enthalpy based equation:

C.4 Mass and momentum conservation model

Non-linear residual

The following residuals represent the weak form of the dilatational flow momentum and mass equations of the arc-plasma.

The following residuals represent the weak form of the incompressible flow momentum and mass equations of the weld pool.

The following residuals represent the weak form of the deformable fluid interface equations of the arc-pool interface.

where l is the element number at the interface, and L is the total number of elements at the interface.

Tangent matrix

The following tangent matrices represent those of the dilatational flow momentum and mass numerical model of the arc-plasma. [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF])

The following tangent matrices represent those of the incompressible flow momentum and mass numerical model of the weld pool. [START_REF] Cuvelier | Some numerical methods for the computation of capillary free boundaries governed by the navier-stokes equations[END_REF])

The following matrices represent those of the deformable fluid interface of the arc-pool interface, projected along the vertical ẑ axis.

Appendix D

Scientific contributions

Peer reviewed articles:

Oral conferences and presentations:

1. 23/09/2020 -Numerical modelling and simulation of welding and additive manufacturing.

Representative of SEMT/LTA on the numerical simulation of welding.

Conseil Scientifique de la Direction des Energies au CEA-Saclay.