# A new strategy for the numerical modeling of a weld pool

Yassine Saadlaoui, Eric Feulvarch, Alexandre Delache, Jean-Baptiste Leblond, Jean-Michel Bergheau

<figure><img src="image_1.png" /></figure>

## To cite this version:

Yassine Saadlaoui, Eric Feulvarch, Alexandre Delache, Jean-Baptiste Leblond, Jean-Michel Bergheau. A new strategy for the numerical modeling of a weld pool. Comptes Rendus Mécanique, 2018, 346 (11), pp.999-1017. 10.1016/j.crme.2018.08.007. hal-01902116

HAL Id: hal-01902116

https://hal.sorbonne-universite.fr/hal-01902116

Submitted on 23 Oct 2018

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Computational methods in welding and additive manufacturing/Simulation numérique des procédés de soudage et de fabrication additive

# A new strategy for the numerical modeling of a weld pool

<figure><img src="image_3.png" /></figure>

Check for
updates

Yassine Saadlaoui<sup>a</sup>, Éric Feulvarch<sup>b,*</sup>, Alexandre Delache<sup>c</sup>,
Jean-Baptiste Leblond<sup>d</sup>, Jean-Michel Bergheau<sup>b</sup>

## ARTICLE INFO

**Article history:**

Received 6 March 2018

Accepted 15 May 2018

Available online 23 August 2018

**Keywords:**

Surface tension

Free surface

Finite-element method

Welding process

## ABSTRACT

Welding processes involve high temperatures and metallurgical and mechanical consequences that must be controlled. For this purpose, numerical simulations have been developed to study the effects of the process on the final structure. During the welding process, the material undergoes thermal cycles that can generate different physical phenomena, like phase changes, microstructure changes and residual stresses and distortions. But the accurate simulation of transient temperature distributions in the part needs to carefully take account of the fluid flow in the weld pool. The aim of this paper is thus to propose a new approach for such a simulation taking account of surface tension effects (including both the “curvature effect” and the “Marangoni effect”), buoyancy forces and free surface motion.

The proposed approach is validated by two numerical tests from the literature: a sloshing test and a plate subjected to a static heat source. Then, the effects of the fluid flow on temperature distributions are discussed in a hybrid laser/arc welding example.

© 2018 Académie des sciences. Published by Elsevier Masson SAS. This is an open access article under the CC BY-NC-ND license ([http://creativecommons.org/licenses/by-nc-nd/4.0/](http://creativecommons.org/licenses/by-nc-nd/4.0/)).

## 1. Introduction

Welding is a very old assembly process dating from the metal age. Different techniques have been developed over time: assemblies by heating and hammering at the beginning, arc welding, resistance welding and more recently, electron beam welding, and laser welding. The large use of welding processes and the need to improve them to achieve a better control of the welding quality led to the development of several numerical simulation methods, which may provide valuable assistance to study the effects of welding on the final structure.

Depending on the objectives targeted by the simulation, different approaches have been developed. These approaches can be classified into two groups, comprising respectively thermo-mechanical simulations to predict the metallurgical and mechanical consequences induced by welding, and simulations of the process itself including the fluid flow in the weld pool. Industrial companies took a large part in the development of thermo-mechanical simulations since the objectives of these simulations are to estimate the residual stresses and distortions [1–5], thus providing very useful information for fatigue lifetime predictions or welding sequence optimization [6]. The thermal part of this kind of simulation is generally simplified

<footnote><sup>a</sup> Université de Lyon, ENISE, UJM, LTDS, UMR 5513 CNRS, 58, rue Jean-Parot, 42023 Saint-Étienne cedex 2, France</footnote>

<footnote><sup>b</sup> Université de Lyon, ENISE, LTDS, UMR 5513 CNRS, 58, rue Jean-Parot, 42023 Saint-Étienne cedex 2, France</footnote>

<footnote><sup>c</sup> LMFA UMR 5509 CNRS, site de Saint-Étienne, Université de Lyon, Université Jean-Monnet de Saint-Étienne, France</footnote>

<footnote><sup>d</sup> UPMC, Université Paris-6, UMR 7190, Institut Jean-Le-Rond-d'Alembert, Tour 65-55, 4, place Jussieu, 75252 Paris cedex 05, France</footnote>

<footnote>* Corresponding author.
E-mail address: eric.feulvarch@enise.fr (É. Feulvarch).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** (a) Simplified schematic of the welding principle, (b) simplified schematic of the material states and fluid flow in the molten pool during the welding process [23].</figcaption></figure>

in order to reduce the calculation time since it is used to simulate large structures [7–9]. The heat input is then represented by an equivalent heat source, the parameters of which must be previously identified on experimental heat affected zones. The thermo-mechanical computation is therefore achieved in two steps: a coupled thermal and metallurgical simulation, followed by a mechanical computation taking account of the temperatures and metallurgical phases previously computed [10,11].

However, this type of simulation is not totally predictive as far as temperature distributions are concerned. Their prediction indeed requires carefully simulating the fluid flow, which plays an important role, particularly on the morphology of the weld pool.

Several authors have therefore developed numerical approaches [12–22], the aim of which is to simulate the physical phenomena involved by the process itself (Fig. 1). Among these phenomena, the fluid flow in the molten pool, particularly related to the surface tension, has a direct effect on the field temperature and weld pool morphology [24–26]. This is why several studies were carried out to model the surface tension in the molten pool [27,28]. The effect of the surface tension on the fluid flow arises from two factors: first, surface tension generates a normal force onto the free surface of the fluid, proportional to the curvature of this surface (“curvature effect”), and second it generates, through a temperature gradient, an additional tangential force onto the free surface (“Marangoni effect”). The classical implementation methods of the surface tension in finite-element codes consider these two effects in a natural way. They include the loads applied on the free surface of the fluid among the external forces of the problem [29,30]. In this type of approach, the implementation of the tangential effect does not raise any problem. However, the normal effect is very difficult to incorporate, since it requires a direct calculation of this force that depends on the free surface motion [31,32].

Several methods have been developed in order to model the free surface during the finite-element computation. These methods can be classified into two categories: the Lagrangian approaches (Arbitrary Lagrangian Eulerian: ALE) [33] and Eulerian approaches (volume of fluid: VOF [34] and level-set [35]). The two Eulerian methods can compute the “keyhole collapse” and help predict defect formation such as porosity [36,37]. However, taking into account the solid deformation using these two methods is very difficult. For this, the ALE approach is often used during the numerical simulation of welding [38–40]. It consists in calculating the displacement of the boundary nodes physically and that of the bulk nodes arbitrarily. This means that the solutions inside the domain are computed with an Eulerian formulation. Thus, the mesh movement is modified in order to improve the mesh quality.

The aim of this paper is to propose a new approach for the simulation of the weld pool taking account of the surface tension effects (including both the “curvature effect” and the “Marangoni effect”), buoyancy forces, and free surface motion.
This simulation will be used in further studies to simulate the interaction between the fluid flow and the solid deformation. For this purpose, the approach developed by Leblond et al. [41] to model the surface tension is used. It has been applied up to now to simulate simple academic benchmarks. In this work, a first application to an industrial problem is proposed (welding process). This modeling of surface tension consists in incorporating the normal and tangential effects in a very simple and efficient way compared to classical methods. In addition, a pragmatic ALE approach is developed to simulate the free surface motion. The choice of this approach is based on two criteria: the simplicity of its implementation, and the possibility to take into account the solid deformation in a simple way.

The paper is organized as follows:

- in section 2, the modeling of the thermal-fluid problem is presented, it contains three parts: an illustration of the fluid and thermal formulation, the modeling of the surface tension, and the modeling of the free surface;

- in section 3, a validation of the modeling of the surface tension and free surface is carried out using numerical validation tests;

- in section 4, the results of the thermal-fluid simulation of laser/arc hybrid welding are presented.

## 2. Modeling the thermal-fluid problem

### 2.1. Strong formulation

Let us consider a bounded fluid domain $\Omega$ of boundary $\partial\Omega$. The coupled heat transfer and fluid flow problems are governed by the following three equations.

**Mass conservation equation**

The fluid is supposed thermally dilatable and mechanically incompressible: the volume variation is only due to temperature variations, the stresses have no effect. From the Eulerian strain rate $D = \frac{1}{2}(\text{grad}\mathbf{v} + (\text{grad}\mathbf{v})^T)$, let us define the Eulerian mechanical strain rate by the formulation:

$$D_m = D + \frac{\dot{\rho}}{3\rho}I \qquad (1)$$

The mass conservation equation then writes:

$$\mathrm{div}\mathbf{v} + \frac{\dot{\rho}}{\rho} = \mathrm{tr}(D_m) = 0 \qquad (2)$$

In practice, the variations of $\rho$ are only due to the thermal dilation effects. It follows that $\frac{\dot{\rho}}{\rho} \approx -\beta\dot{T}$, where $\beta$ is the coefficient of volumetric expansion.

**Energy balance equation**

$$\rho \left( \frac{\partial H}{\partial t} + \mathbf{v} \cdot \text{grad} H \right) - \mathrm{div}(\lambda \text{grad} T) - r = 0 \qquad (3)$$

Here $\rho$ is the density, $H$, the specific enthalpy, $\mathbf{v}$, the fluid velocity, $\lambda$, the heat conductivity, $T$, the temperature, and $r$ the volume heat source. The enthalpy $H$ is a function of temperature including both the sensible heat and the latent heat due to phase change, and $\rho$ and $\lambda$ can be temperature-dependent. The associated boundary conditions are:

$$T = T_0 \quad \text{on} \quad \partial\Omega_T \qquad (4)$$

$$\lambda \text{grad} T \cdot \mathbf{n} = q \quad \text{on} \quad \partial\Omega_q \qquad (5)$$

$T_0$ is the known temperature on $\partial\Omega_T$ and $q$ a temperature-dependent heat flux on $\partial\Omega_q$ whose outside unit normal is $\mathbf{n}$.

**Momentum equations (Navier-Stokes equations)**

$$-\text{grad} p + \mathrm{div}\mathbf{s} + \mathbf{f}_v = \rho \left( \frac{\partial \mathbf{v}}{\partial t} + \text{grad} \mathbf{v} \cdot \mathbf{v} \right) \qquad (6)$$

The fluid flow is supposed to be laminar. Here $p$ is the fluid pressure, $\mathbf{s}$, the stress deviator, and $\mathbf{f}_v$, the volume force. Buoyancy forces can be considered using the volume forces which are given by the following equation:

$$\mathbf{f}_v = \rho_0(1 - \beta(T - T_0))\mathbf{g} \qquad (7)$$

$\mathbf{g}$ being the gravitational acceleration and $\rho_0$ the mass per unit volume at the reference temperature $T_0$.

The fluid is supposed to behave as a Newtonian fluid. The relation between the Eulerian mechanical strain rate $D_m$ and the stress deviator $\mathbf{s}$ can therefore be written as:
<figure><img src="image_5.png" /><figcaption>Fig. 2. P1 + /P1 element.</figcaption></figure>

$$ \mathbf{s} = 2\mu\mathbf{D}_m \qquad (8) $$

where $\mu$ is the fluid viscosity.

The associated boundary conditions are:

$$ \mathbf{v} = \mathbf{v}^d \quad \text{on } \partial\Omega_v \qquad (9) $$

$$ \sigma \cdot \mathbf{n} = \varphi^d \quad \text{on } \partial\Omega_\varphi \qquad (10) $$

Here $\mathbf{v}^d$ and $\varphi^d$ are the prescribed velocity and surface force on $\partial\Omega_v$ and $\partial\Omega_\varphi$ respectively, and $\sigma = \mathbf{s} - p\mathbf{I}$ is the Cauchy stress tensor.

## 2.2. Finite element formulation

The weak formulation of the energy balance equation as well as the finite element formulation of the thermal problem are quite classical, and details can be found in [42,43]. We therefore focus on the fluid problem formulation in the sequel.

The finite element formulation of the fluid problem rests on a mixed velocity-pressure ($\mathbf{v}, p$) variational formulation.

The problem is now to find $(\mathbf{v}, p) \in (V_{ad} \times P_{ad}), \forall t \in [0, T]$, such as $\forall (\mathbf{v}^*, p^*) \in (V_{ad}^0 \times P_{ad})$:

$$ - \int_{\Omega} \mathbf{v}^* \rho \left( \frac{\partial \mathbf{v}}{\partial t} + \mathbf{grad} \mathbf{v} \cdot \mathbf{v} \right) d\mathbf{v} - \int_{\Omega} \mathbf{D}^* : 2\mu \mathbf{D}_m d\mathbf{v} + \int_{\Omega} \mathrm{div}(\mathbf{v}^*) p d\mathbf{v} + \int_{\Omega} \mathbf{v}^* \cdot \mathbf{f}_v d\mathbf{v} + \int_{\partial\Omega_\varphi} \mathbf{v}^* \cdot \varphi^d d\mathbf{s} = 0 \qquad (11) $$

$$ \int_{\Omega} p^* \left( \mathrm{div} \mathbf{v} + \frac{\dot{\rho}}{\rho} \right) d\mathbf{v} = 0 $$

with the initial condition $\mathbf{v}(\mathbf{x}, t=0) = \mathbf{v}_0(\mathbf{x})$.

In Eqs. (11), $\mathbf{D}^* = \frac{1}{2}(\mathbf{grad} \mathbf{v}^* + (\mathbf{grad} \mathbf{v}^*)^T)$ and:

$$ V_{ad} = \{ \mathbf{v} \in H^1(\Omega) \mid \mathbf{v} = \mathbf{v}^d \text{ on } \partial\Omega_v \} $$

$$ V_{ad}^0 = \{ \mathbf{v} \in H^1(\Omega) \mid \mathbf{v} = \mathbf{0} \text{ on } \partial\Omega_v \} $$

$$ P_{ad} = \{ p \in L^2(\Omega) \} $$

P1 + /P1 tetrahedral elements (Fig. 2) are used for the discretization of Eqs. (11).

### Velocity approximation

The velocity field is decomposed into a part $\mathbf{v}_s$ calculated from the velocities of the vertex nodes gathered in the array $\mathbf{v}_s$, using the shape functions $\mathbf{N}_s$ of the classical four-node tetrahedron [42], plus another part $\mathbf{v}_b$ calculated from a vector $\mathbf{\Lambda}_b$ of additional degrees of freedom of the bubble (internal) node and the associated shape function $\mathbf{N}_b$:

$$ \mathbf{v} = \mathbf{v}_s + \mathbf{v}_b \qquad (12) $$

with

$$ \mathbf{v}_s = \mathbf{N}_s \cdot \mathbf{v}_s = \sum_{n=1}^{4} \mathbf{N}_{s_n} \mathbf{v}_{s_n} \qquad (13) $$

$$ \mathbf{v}_b = \mathbf{N}_b \cdot \mathbf{\Lambda}_b $$

with $\mathbf{N}_b = \min_{n=1}^4 \mathbf{N}_{s_n}$.

Note that $\mathbf{\Lambda}_b$ does not represent the velocity of the centroid of the element because the shape functions of the vertex nodes $\mathbf{N}_s$ do not vanish at this point.
The virtual velocity $\mathbf{v}^*$ is decomposed in the same way. Note that the “bubble velocity” $\mathbf{v}_b$ vanishes on each element boundary $\partial\Omega^e$ (bubble node principle) and therefore at any point of $\partial\Omega$.

### Pressure approximation

The pressure is calculated from the values of the vertex nodes $\mathbf{P}_s$:

$$p = \mathbf{N}_s \cdot \mathbf{P}_s = \sum_{\mathbf{n}=1}^{4} N_{s\mathbf{n}} P_{s\mathbf{n}} \qquad (14)$$

and so is also the temperature.

The Eulerian mechanical strain rate $\mathbf{D}_m$, as well as its virtual counterpart $\mathbf{D}_m^*$, are then decomposed in the following form:

$$\begin{align}
\mathbf{D}_m &= \mathbf{D}_{m_s} + \mathbf{D}_{m_b} \nonumber \\
\mathbf{D}_{m_s} &= \mathbf{D}_s + \frac{\dot{\rho}}{3\rho} \mathbf{I} = \frac{1}{2} (\mathbf{grad} \mathbf{v}_s + (\mathbf{grad} \mathbf{v}_s)^\mathrm{T}) + \frac{\dot{\rho}}{3\rho} \mathbf{I} \tag{15} \\
\mathbf{D}_{m_b} &= \mathbf{D}_b = \frac{1}{2} (\mathbf{grad} \mathbf{v}_b + (\mathbf{grad} \mathbf{v}_b)^\mathrm{T}) \nonumber
\end{align}$$

Equations (11) with (8), (12) and (15) can therefore be rewritten, for any $\mathbf{v}_s^*, \mathbf{v}_b^*$, and $p^*$:

$$\begin{align}
& -\int_{\Omega} \mathbf{v}_s^* \rho \left( \frac{\partial \mathbf{v}}{\partial t} + \mathbf{grad} \mathbf{v} \cdot \mathbf{v} \right) d\mathbf{v} - \int_{\Omega} \mathbf{D}_s^* : 2\mu \mathbf{D}_{m_s} d\mathbf{v} - \int_{\Omega} \mathbf{D}_s^* : 2\mu \mathbf{D}_{m_b} d\mathbf{v} + \int_{\Omega} \mathrm{div}(\mathbf{v}_s^*) p d\mathbf{v} \nonumber \\
& \qquad + \int_{\Omega} \mathbf{v}_s^* \cdot \mathbf{f}_v d\mathbf{v} + \int_{\partial\Omega_\varphi} \mathbf{v}_s^* \cdot \boldsymbol{\varphi}^\mathrm{d} d\mathbf{s} = 0 \nonumber \\
& -\int_{\Omega} \mathbf{v}_b^* \rho \left( \frac{\partial \mathbf{v}}{\partial t} + \mathbf{grad} \mathbf{v} \cdot \mathbf{v} \right) d\mathbf{v} - \int_{\Omega} \mathbf{D}_b^* : 2\mu \mathbf{D}_{m_s} d\mathbf{v} - \int_{\Omega} \mathbf{D}_b^* : 2\mu \mathbf{D}_{m_b} d\mathbf{v} + \int_{\Omega} \mathrm{div}(\mathbf{v}_b^*) p d\mathbf{v} \tag{16} \\
& \qquad + \int_{\Omega} \mathbf{v}_b^* \cdot \mathbf{f}_v d\mathbf{v} = 0 \nonumber \\
& \int_{\Omega} p^* \left( \mathrm{div}(\mathbf{v}_s + \mathbf{v}_b) + \frac{\dot{\rho}}{\rho} \right) d\mathbf{v} = 0 \nonumber
\end{align}$$

$\mathbf{D}_s$ and $\mathbf{D}_s^*$ are constant over each element $\Omega^e$ as they involve the derivatives of the shape functions of the classical four-node tetrahedron. Assuming that $\frac{\dot{\rho}}{\rho}$ and $\mu$ are constant over each element (in practice $\frac{\dot{\rho}}{\rho}$ and $\mu$ are calculated at the centroid of the tetrahedron), $\mathbf{D}_{m_s}$ is also constant over each element. Since $\mathbf{v}_b$ vanishes on the element boundary $\partial\Omega^e$, it follows that:

$$\begin{align}
\int_{\Omega^e} \mathbf{D}_s^* : 2\mu \mathbf{D}_{m_b} d\mathbf{v} &= 2\mu \mathbf{D}_s^* : \int_{\Omega^e} \mathbf{D}_{m_b} d\mathbf{v} = \mu \mathbf{D}_s^* : \int_{\Omega^e} (\mathbf{grad} \mathbf{v}_b + (\mathbf{grad} \mathbf{v}_b)^\mathrm{T}) d\mathbf{v} \nonumber \\
&= \mu \mathbf{D}_s^* : \int_{\partial\Omega^e} (\mathbf{n} \otimes \mathbf{v}_b + \mathbf{v}_b \otimes \mathbf{n}) d\mathbf{s} = 0 \tag{17}
\end{align}$$

In the same way, one can show that:

$$\int_{\Omega^e} \mathbf{D}_b^* : 2\mu \mathbf{D}_{m_s} d\mathbf{v} = 0 \qquad (18)$$

An implicit backward scheme for $\frac{\partial \mathbf{v}}{\partial t}$ and a linearization of the convection term are used for the time integration of Eqs. (16). It is also supposed that the bubble node does not contribute to the acceleration (just the vertex nodes are taken into account for the interpolation of acceleration).

Using Eqs. (17) and (18), Eqs. (16) can therefore be rewritten at the final instant of any time-step $[t, t + \Delta t]$, for any $\mathbf{v}_s^*, \mathbf{v}_b^*$, and $p^*$:
$$
\begin{equation}
\begin{aligned}
& - \int_{\Omega} \mathbf{v}_s^* \rho \left( \frac{\mathbf{v}_s - \mathbf{v}_s^t}{\Delta t} + \mathbf{grad}(\mathbf{v}_s + \mathbf{v}_b) \cdot \mathbf{v}^t \right) d\mathbf{v} - \int_{\Omega} \mathbf{D}_s^* : 2\mu \mathbf{D}_{ms} d\mathbf{v} + \int_{\Omega} \mathrm{div}(\mathbf{v}_s^*) p d\mathbf{v} \\
& \qquad + \int_{\Omega} \mathbf{v}_s^* \cdot \mathbf{f}_v d\mathbf{v} + \int_{\partial \Omega_{\varphi}} \mathbf{v}_s^* \cdot \varphi^d d\mathbf{s} = 0 \\
& - \int_{\Omega} \mathbf{v}_b^* \rho \left( \frac{\mathbf{v}_b - \mathbf{v}_b^t}{\Delta t} + \mathbf{grad}(\mathbf{v}_b + \mathbf{v}_b) \cdot \mathbf{v}^t \right) d\mathbf{v} - \int_{\Omega} \mathbf{D}_b^* : 2\mu \mathbf{D}_{mb} d\mathbf{v} + \int_{\Omega} \mathrm{div}(\mathbf{v}_b^*) p d\mathbf{v} + \int_{\Omega} \mathbf{v}_b^* \cdot \mathbf{f}_v d\mathbf{v} = 0 \\
& \int_{\Omega} p^* \left( \mathrm{div}(\mathbf{v}_s + \mathbf{v}_b) + \frac{\dot{\rho}}{\rho} \right) d\mathbf{v} = 0
\end{aligned}
\tag{19}
\end{equation}
$$

where all the quantities are considered at time $t + \Delta t$. except for $\mathbf{v}^t$, which represents the velocity at time $t$.

Using now the Voigt notation for the Eulerian strain rate and defining, from (13) and (15), $\mathbf{B}_s$ and $\mathbf{B}_b$ matrices such that $\mathbf{D}_s = \mathbf{B}_s \cdot \mathbf{V}_s$ and $\mathbf{D}_b = \mathbf{B}_b \cdot \mathbf{\Lambda}_b$, Eqs. (19) lead to the following expressions of the nodal residual vectors:

$$
\begin{align}
\mathbf{T}_s &= -\int_{\Omega} \mathbf{N}_s^\mathrm{T} \cdot \frac{\mathbf{v}_s - \mathbf{v}_s^t}{\Delta t} \rho \, d\mathbf{v} - \int_{\Omega} \mathbf{N}_s^\mathrm{T} \cdot \rho \, \mathbf{grad}(\mathbf{v}_s + \mathbf{v}_b) \cdot \mathbf{v}^t \, d\mathbf{v} \nonumber \\
&\quad - \int_{\Omega} \mathbf{B}_s^\mathrm{T} \cdot 2\mu \mathbf{D}_{ms} \, d\mathbf{v} + \int_{\partial \Omega_\varphi} (\mathbf{grad} \, \mathbf{N}_s)^\mathrm{T} \, p \, d\mathbf{v} + \int_{\Omega} \mathbf{N}_s^\mathrm{T} \cdot \varphi^d \, d\mathbf{s} + \int_{\Omega} \mathbf{N}_s^\mathrm{T} \cdot \mathbf{f}_v \, d\mathbf{v} = 0 \tag{20} \\
\mathbf{T}_b &= -\int_{\Omega} \mathbf{N}_b^\mathrm{T} \cdot \frac{\mathbf{v}_b - \mathbf{v}_b^t}{\Delta t} \rho \, d\mathbf{v} - \int_{\Omega} \mathbf{N}_b^\mathrm{T} \cdot \rho \, \mathbf{grad}(\mathbf{v}_s + \mathbf{v}_b) \cdot \mathbf{v}^t \, d\mathbf{v} \nonumber \\
&\quad - \int_{\Omega} \mathbf{B}_b^\mathrm{T} \cdot 2\mu \mathbf{D}_{mb} \, d\mathbf{v} + \int_{\partial \Omega_\varphi} (\mathbf{grad} \, \mathbf{N}_b)^\mathrm{T} \, p \, d\mathbf{v} + \int_{\Omega} \mathbf{N}_b^\mathrm{T} \cdot \mathbf{f}_v \, d\mathbf{v} = 0 \nonumber \\
\mathbf{T}_p &= \int_{\Omega} \mathbf{N}_s^\mathrm{T} \left( \mathrm{div}(\mathbf{v}_s + \mathbf{v}_b) + \frac{\dot{\rho}}{\rho} \right) \, d\mathbf{v} = 0 \nonumber
\end{align}
$$

Note that the second equation of system (20) defines, at the element level, the following linear relation between $\mathbf{V}_s$, $\mathbf{\Lambda}_b$ and $\mathbf{P}_s$:

$$
\mathbf{K}_{bs} \cdot \mathbf{V}_s + \mathbf{K}_{bb} \cdot \mathbf{\Lambda}_b + \mathbf{K}_{bp} \cdot \mathbf{P}_s - \mathbf{F}_b = 0 \quad (21)
$$

with:

$$
\begin{align}
\mathbf{K}_{bs} &= \int_{\Omega} \mathbf{N}_{b}^{\mathrm{T}} \cdot \rho \mathbf{grad} \, \mathbf{N}_{s} \cdot \mathbf{v}^{t} d\mathbf{v} + \int_{\Omega} \mathbf{N}_{b}^{\mathrm{T}} \cdot \frac{\rho}{\Delta t} \mathbf{N}_{s} d\mathbf{v} \nonumber \\
\mathbf{K}_{bb} &= \int_{\Omega} \mathbf{B}_{b}^{\mathrm{T}} \cdot 2\mu \mathbf{B}_{b} d\mathbf{v} + \int_{\Omega} \mathbf{N}_{b}^{\mathrm{T}} \cdot \rho \mathbf{grad} \, \mathbf{N}_{b} \cdot \mathbf{v}^{t} d\mathbf{v} \nonumber \\
\mathbf{K}_{bp} &= - \int_{\Omega} (\mathbf{grad} \, \mathbf{N}_{b})^{\mathrm{T}} \cdot \mathbf{N}_{s} d\mathbf{x} \tag{22} \\
\mathbf{F}_{b} &= \int_{\Omega} \mathbf{N}_{b}^{\mathrm{T}} \cdot \mathbf{f}_{v} d\mathbf{v} + \int_{\Omega} \mathbf{N}_{b}^{\mathrm{T}} \cdot \frac{\mathbf{v}_{s}^{t}}{\Delta t} \rho d\mathbf{v} \nonumber
\end{align}
$$

Equations (20) are solved using a Newton-Raphson method. After elimination of $\mathbf{\Lambda}_b$ from Eq. (21), at each iteration corrections $\delta \mathbf{V}_s$ and $\delta \mathbf{P}_s$ are calculated through the following equations:

$$
\begin{bmatrix} \mathbf{K}_{ss} - \mathbf{K}_{sb} \cdot \mathbf{K}_{bb}^{-1} \cdot \mathbf{K}_{bs} & \mathbf{K}_{sp} - \mathbf{K}_{sb} \cdot \mathbf{K}_{bb}^{-1} \cdot \mathbf{K}_{bp} \\ \mathbf{K}_{ps} - \mathbf{K}_{pb} \cdot \mathbf{K}_{bb}^{-1} \cdot \mathbf{K}_{bs} & -\mathbf{K}_{pb} \cdot \mathbf{K}_{bb}^{-1} \cdot \mathbf{K}_{bp} \end{bmatrix} \begin{Bmatrix} \delta \mathbf{V}_s \\ \delta \mathbf{P}_s \end{Bmatrix} = \begin{Bmatrix} \mathbf{T}_s \\ \mathbf{T}_p \end{Bmatrix} \quad (23)
$$

with:
$$
\begin{align}
\mathbf{K}_{sb} &= \int_{\Omega} \mathbf{N}_{\mathbf{s}}^{\mathrm{T}} \cdot \rho \, \mathrm{grad} \, \mathbf{N}_{\mathbf{b}} \cdot \mathbf{v}^{\mathrm{T}} \, \mathrm{d}v \nonumber \\
\mathbf{K}_{ss} &= \int_{\Omega} \mathbf{B}_{\mathbf{s}}^{\mathrm{T}} \cdot 2\mu \mathbf{B}_{\mathbf{s}} \, \mathrm{d}v + \int_{\Omega} \mathbf{N}_{\mathbf{s}}^{\mathrm{T}} \cdot \rho \, \mathrm{grad} \, \mathbf{N}_{\mathbf{s}} \cdot \mathbf{v}^{\mathrm{T}} \, \mathrm{d}v + \int_{\Omega} \mathbf{N}_{\mathbf{s}}^{\mathrm{T}} \cdot \frac{\rho}{\Delta t} \, \mathrm{grad} \, \mathbf{N}_{\mathbf{s}} \, \mathrm{d}v \tag{24} \\
\mathbf{K}_{sp} &= - \int_{\Omega} (\mathrm{grad} \, \mathbf{N}_{\mathbf{s}})^{\mathrm{T}} p \cdot \mathbf{N}_{\mathbf{s}} \, \mathrm{d}x \nonumber \\
\mathbf{K}_{ps} &= \mathbf{K}_{sp}^{\mathrm{T}} \nonumber \\
\mathbf{K}_{pb} &= \mathbf{K}_{bp}^{\mathrm{T}} \tag{25}
\end{align}
$$

The process is repeated as long as $||\mathbf{T}_s|| > \varepsilon_s$ and $||\mathbf{T}_p|| > \varepsilon_p$, where $\varepsilon_s$ and $\varepsilon_p$ are precision thresholds prescribed by the user.

## 2.3. Modeling the surface tension effects

One of the main effects to be accounted for in weld pool simulations is related to surface tension. Surface tension induces a pressure normal to the free surface of the fluid, proportional to the curvature of this surface, and its gradient generates a tangential surface force (Marangoni effect). As far as the implementation of the surface tension in finite-element codes [44,45] is concerned, these two effects are taken into account in a natural way through the boundary conditions as defined in Eq. (7) by including the loads applied on the free surface of the fluid among the external forces of the problem. In such an approach, if the Marangoni effect can be introduced in a straightforward way, the implementation of the curvature effect is far more complicated, as it requires the preliminary calculation of the mean curvature (trace of the curvature tensor) of the free surface [31,32]. But the position of the free surface is given by the fluid motion, and the curvature calculation requires the evaluation of the spatial second derivatives of the displacement field. This calculation can be achieved quite accurately using $C^1$ elements, but these elements are very heavy and time-consuming. When classical $C^0$ elements are used, another way consists in calculating the divergence of the unit normal to the free surface, and the spatial first derivatives of the shape functions [46,47]. But the calculation of the unit normal, generally obtained by averaging the normal of the surface elements connected to each node considered, can induce significant errors in the results.

In this study, we use the efficient approach suggested in [27,44,45] for the curvature effect and generalized in [41] for both the curvature and Marangoni effects. The approach is based on the equivalence between the effects of the surface tension and those of a fictitious stressed membrane bonded to the free surface. The surface tension is therefore taken into account through a surface hydrostatic state of stress $\boldsymbol{\tau} = \begin{bmatrix} \gamma & 0 \\ 0 & \gamma \end{bmatrix}$, where $\gamma$ is the surface tension, applied to membrane-type (skin) elements (Fig. 3) that do not have any material properties [41]. The transmission of the surface tension effects from the skin elements to the 3D fluid elements is carried out naturally through the common nodes of the 2D and 3D elements.

The weak formulation of the fluid problem is now modified by adding the term $-\int_{\partial\Omega_{\text{free}}} \mathbf{D}_{\mathbf{s}}^* : \boldsymbol{\tau} \, \mathrm{d}s$, where $\partial\Omega_{\text{free}}$ is the free surface of the fluid, in the left-hand side of the first of Eq. (19), which is now written:

$$
\begin{equation}
\begin{split}
& -\int_{\Omega} \mathbf{v}_{\mathbf{s}}^{*} \rho \left( \frac{\mathbf{v}_{\mathbf{s}} - \mathbf{v}_{\mathbf{s}}^{\mathrm{T}}}{\Delta t} + \mathbf{grad}(\mathbf{v}_{\mathbf{s}} + \mathbf{v}_{\mathbf{b}}) \cdot \mathbf{v}^{\mathrm{T}} \right) \mathrm{d}v - \int_{\Omega} \mathbf{D}_{\mathbf{s}}^{*} : 2\mu \mathbf{D}_{\mathbf{m}_{\mathbf{s}}} \mathrm{d}v + \int_{\Omega} \mathrm{div}(\mathbf{v}_{\mathbf{s}}^{*}) p \mathrm{d}v \\
& - \int_{\partial\Omega_{\text{free}}} \mathbf{D}_{\mathbf{s}}^{*} : \boldsymbol{\tau} \mathrm{d}s + \int_{\Omega} \mathbf{v}_{\mathbf{s}}^{*} \cdot \mathbf{f}_{\mathbf{v}} \mathrm{d}v + \int_{\partial\Omega_{\varphi}} \mathbf{v}_{\mathbf{s}}^{*} \cdot \varphi^{\mathrm{d}} \mathrm{d}s = 0
\end{split}
\tag{26}
\end{equation}
$$

The proposed approach does not require calculating the mean curvature of the free surface, which makes it efficient and accurate. In addition, it is possible to account for both normal and tangential effects in a natural and simple way by simply considering temperature-dependent values of $\gamma$. It is sufficient to consider the correct value of the surface tension at each Gauss point of the membrane elements. This allows us to avoid the numerical problems of classical approaches.

## 2.4. Modeling the free surface evolution

To take the free surface evolution into account, different methods have been proposed. Combined with fully Eulerian formulations, the volume of fluid [34] or level set [35] techniques are generally used. Alternatively, the Arbitrary Lagrangian–Eulerian (ALE) approach consists in imposing the free surface motion to the mesh [38–40]. But in order to keep a good mesh quality, the nodes inside the weld pool must in turn be repositioned, thus defining a mesh motion which must be taken into account in the finite element formulation.

The algorithm used in the present work is summarized in Fig. 4. The first step consists in recovering the mesh and the results for a given time step of the thermal-fluid computation. Here, a strong coupling between fluid and thermal
<figure><img src="image_6.png" /><figcaption>Fig. 3. Mesh containing skin elements (red 1D elements with sky blue nodes).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Proposed algorithm of the ALE approach.</figcaption></figure>

computations is used. Then the status of each node (free surface, fluid, solid) is defined. The free-surface nodes are then displaced in order to follow the free surface motion. The fluid nodes inside the molten pool are finally repositioned using an elastic calculation where the solid nodes remain fixed.

The elements having a mean temperature above the melting temperature $T_{liq}$ are considered as fluid elements, the others as solid elements. For each node, depending on its position and the status of the elements it is connected to, four-node statuses can be defined:

- solid node (red node): connected to at least one solid element;

- solid/free-surface node (purple node) – solid node connected to one of the skin elements defining the free surface;

- fluid node (green node): connected only to fluid elements;

- fluid/free-surface node (sky blue node): fluid node connected to one of the skin elements defining the free surface.

The normal at a free-surface node is defined as the mean value of the normal of the skin elements containing it (Fig. 5).

The normal velocities and the displacements of the free-surface nodes are then calculated using the following expressions:

$$ \mathbf{v}_{n_k}(t + \Delta t) = (\mathbf{n}_k \cdot \mathbf{v}_k(t + \Delta t)) \cdot \mathbf{n}_k \qquad (27) $$
<figure><img src="image_8.png" /><figcaption>Fig. 5. Definition of the normal at a node of the free surface and projection of material velocity.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. The mesh before updating the node positions.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. The mesh after updating the node positions.</figcaption></figure>

$$ \mathbf{u}_{n_k}(t + \Delta t) = \mathbf{u}_{n_k}(t) + \mathbf{v}_{n_k}(t + \Delta t) \cdot \Delta t \quad (28) $$

Here $\mathbf{v}_{n_k}$ and $\mathbf{u}_{n_k}$ are respectively the normal velocity and the normal displacement of the free-surface node $k$, $\mathbf{v}_k$ is the velocity of this node, and $\mathbf{n}_k$ is its normal (Fig. 6).

Once the displacements of the free-surface nodes have been calculated, different methods can be used to reposition the internal fluid nodes. For example, Heuzé et al. [48] used a barycentric replacement. In the present approach, an elastic computation is carried out in order to find the arbitrary displacements of the fluid nodes. In this computation, the displacements of the solid nodes and the free-surface nodes are prescribed. Finally, the mesh is updated using the computed displacements (Fig. 7).

All these steps are performed at each time step of the thermal-fluid computation.

**Remark.** For the sake of simplicity, the mesh velocity has not been included in the presentation of the finite element formulation above. But of course, the mesh velocity is taken into account in all the advection terms such as $\text{grad}(\mathbf{v}_s + \mathbf{v}_b) \cdot \mathbf{v}^t$ in Eq. (19), which is now written $\text{grad}(\mathbf{v}_s + \mathbf{v}_b) \cdot (\mathbf{v}^t - \mathbf{v}_{\text{mesh}}^t)$, where $\mathbf{v}_{\text{mesh}}^t$ is the mesh velocity.
<figure><img src="image_11.png" /><figcaption>Fig. 8. Schematization and mesh of the sloshing problem.</figcaption></figure>

## 3. Validation tests

### 3.1. Sloshing problem

The sloshing problem was developed by Ramaswamy et al. [49] to validate an ALE formulation. It consists in following the oscillations of a liquid in a container. Sloshing was initially considered to be due only to gravity. This problem was generalized by other authors by accounting for more physical phenomena. For example, Popinet et al. [50] have simulated in 2D the oscillations produced by the surface tension. They have compared their numerical results with the analytic results given by Prosperetti et al. [51]. The problem was considered later by El-Sayed et al. [25] using a 3D mesh.

In this work, the sloshing problem was used to validate the ALE approach of the free surface and the normal effect of the surface tension. The aim of the numerical model is to follow the oscillations of a single liquid in a container that has a square geometry. Initially, the free surface of the liquid has a sinusoidal shape given by the following equation:

$$H(x) = 1 + a_0 \sin(\pi(0.5 - x)) \qquad (29)$$

where $a_0 = 0.01$ is the initial amplitude. In this equation all distances are expressed in m.

With regard to the boundary conditions, the bottom surface is clamped and the normal velocities are such that $v \cdot n = 0$ on the side walls. The free surface is the top one. The mesh (Fig. 8) contains 4223 nodes and 19,736 linear tetrahedral elements. The free surface is meshed with 100 2D triangular “skin” elements (red elements in Fig. 8) to apply the surface tension. The liquid properties are: relative density $\rho = 1 \text{ kg} \cdot \text{m}^{-3}$, kinematic viscosity $\eta = 0.01 \text{ m}^2 \cdot \text{s}^{-1}$, and gravity $g = 1 \text{ m} \cdot \text{s}^{-2}$ (artificial value).

A first simulation is carried out to follow the sloshing of a liquid under the effect of its sole weight (gravity effect). The results of this simulation and the analytical solution [51] are provided by Fig. 9.

Regarding the other simulation, the sloshing of the liquid is generated by its weight as well as the surface tension. Values of the surface tension between 0.001 and $1 \text{ N} \cdot \text{m}^{-1}$ are used. The results show proportionality between the value of the surface tension and the intensity of liquid sloshing (Fig. 10). The comparison between the numerical results and the analytical solution shows a good correlation: the maximal error is less than 0.006%. We note that the value of this error is very similar to those found in the literature. For example, Rabier et al. [32] have found a maximal error equal to 0.005%. However, they meshed with a discretization twice finer than ours. In addition to that, they used a time step equal to 0.0025 s, that is $\frac{1}{8}$ of the value used in our simulation (0.1 s).

### 3.2. Benchmark problem

The Benchmark consists of some numerical simulation of the spot welding “GTWA” process (Fig. 11). It was developed in a 2D axisymmetric option by Girard et al. [52], in a 3D option by Hamide et al. [53] and El-Sayed et al. [25] to validate the tangential effect of the surface tension. In this work, a purely thermal-fluid computation of this problem was carried out on the quarter of a disk (Fig. 11). A viscosity depending on the temperature allows separating the fluid and solid zones. The solid is thus schematized as a fluid with a high viscosity ($\mu = 10^6 \text{ kg} \cdot \text{m}^{-1} \cdot \text{s}^{-1}$). The mesh includes 68,130 linear tetrahedral elements, 762 skin elements and 13,617 nodes. The length of the elements edges is 30 µm in the fluid zone and more than 100 µm in the solid zone.

The material used is steel, supposed to be isotropic, homogeneous, Newtonian, thermally dilatable, and mechanically incompressible. Its properties are given in Table 1.
<figure><img src="image_12.png" /><figcaption>Fig. 9. The displacement of the free surface under the effect of the liquid weight.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. The displacement of the free surface under the effect of the liquid weight and the surface tension; $\gamma = 0.001 \text{ N}\cdot\text{m}^{-1}$ and $\gamma = 1 \text{ N}\cdot\text{m}^{-1}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Schematization [52] and mesh of the benchmark problem.</figcaption></figure>

With regard to the boundary conditions, the node velocities of the free surface and symmetric boundaries are such that $\mathbf{v} \cdot \mathbf{n} = 0$. The heat flux imposed on the upper surface of the disk is the sum of three terms:

$$ \varnothing_{\text{total}} = \varnothing_{\text{source}} - \varnothing_{\text{convection}} - \varnothing_{\text{radiation}} \quad (30) $$

where:

$$ \varnothing_{\text{source}}(x) = \frac{P}{2\pi R_0^2} \exp\left(-\frac{R^2}{2R_0^2}\right) \quad (31) $$
<table><caption>**Table 1**

Characteristics of the material.</caption><tr><td>Solidus temperature</td><td>T<sub>S</sub> = 1696 K</td></tr><tr><td>Liquidus temperature</td><td>T<sub>L</sub> = 1740 K</td></tr><tr><td>Relative density (solid)</td><td>ρ = 7500 kg·m<sup>-3</sup></td></tr><tr><td>Specific heat (solid)</td><td>c = 602 J·K<sup>-1</sup>·kg<sup>-1</sup></td></tr><tr><td>Thermal conductivity (solid)</td><td>λ = 24 W·m<sup>-1</sup>·K<sup>-1</sup></td></tr><tr><td>Relative density (liquid)</td><td>ρ = 6350 kg·m<sup>-3</sup></td></tr><tr><td>Specific heat (liquid)</td><td>c = 695 J·K<sup>-1</sup>·kg<sup>-1</sup></td></tr><tr><td>Thermal conductivity (liquid)</td><td>λ = 20 W·m<sup>-1</sup>·K<sup>-1</sup></td></tr><tr><td>Dynamic viscosity</td><td>μ = 2.5 · 10<sup>-3</sup> kg·m<sup>-1</sup>·s<sup>-1</sup></td></tr><tr><td>Coefficient of expansion</td><td>β = 10<sup>-4</sup> K<sup>-1</sup></td></tr></table>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Temperature evolution with positive derivative of surface tension, $\frac{\partial\gamma}{\partial T} > 0$.</figcaption></figure>

with $R_0 = 3 \cdot 10^{-3}$ m and $P = 790$ W.

$$ \varnothing_{\text{radiation}} = \epsilon \sigma_0 (T^4 - T_0^4) \qquad (32) $$

with $\sigma_0$ the Stefan-Boltzmann constant, $\epsilon$ the emissivity equal to 0.5, and $T_0 = 300$ K.

$$ \varnothing_{\text{convection}} = h(T - T_0) \qquad (33) $$

with $h$ the transfer coefficient equal to $15 \text{ W} \cdot \text{m}^{-2} \cdot \text{K}^{-1}$.

The lateral edges of the disk are adiabatic, and a linear temperature decrease from 800 to 400 K is enforced on its lower side, from the center to the edge (Fig. 11).

During this simulation, only the tangential effect of the surface tension is taken into account, without the normal effect. Two values of the derivative of the surface tension, $\frac{\partial\gamma}{\partial T} = 10^{-4} \text{ N} \cdot \text{m}^{-1} \cdot \text{K}^{-1}$ and $\frac{\partial\gamma}{\partial T} = -10^{-4} \text{ N} \cdot \text{m}^{-1} \cdot \text{K}^{-1}$, are imposed on the skin elements of the upper face of the disk. The first three seconds of the process are simulated, and the time step is 0.1 s during the formation of the weld pool (thermal computation) and $2 \cdot 10^{-4}$ s after the activation of the fluid flow (thermal-fluid computation).

A comparison between the Benchmark and our numerical results is carried out. It consists in examining the following quantities:

- temperature on the radius of the upper face $T(z=0)$.

- radial velocity $v_r(z=0)$.

Figs. 12 and 13 show respectively the evolution of the temperature field and the velocity vector with a positive derivative of the surface tension with respect to temperature. Figs. 14 and 15 illustrate the case of a negative derivative. A good correlation between our results and the benchmark results is found in both cases. With regard to the benchmark results, we use the average value of the results given by five computation codes [25,53]. A maximal error less than 2% is found for the temperature. Regarding the velocity, a maximal error of 10% is found. The magnitudes of the temperature and velocity are also checked on the symmetry axis and here again, a good correlation with previous simulations is found.
<figure><img src="image_16.png" /><figcaption>Fig. 13. Velocity evolution with positive derivative of surface tension, $\frac{\partial\gamma}{\partial T} > 0$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Temperature evolution with negative derivative of surface tension, $\frac{\partial\gamma}{\partial T} < 0$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 15. Velocity evolution with negative derivative of surface tension, $\frac{\partial\gamma}{\partial T} < 0$.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 16. Finite-element model: mesh.</figcaption></figure>

# 4. Numerical simulation of the hybrid welding process

## 4.1. FE model

In this section, a 3D thermal-fluid simulation of hybrid laser/arc welding is presented using the SYSWELD® code [54]. It takes into account the fluid flow through the surface tension and the buoyancy. Both equivalent double ellipsoidal and Gaussian surface heat sources [55] are applied in the Y direction to a metal block 30 mm in length, 10 mm in width, and 20 mm in height. The model has a symmetry about the YZ plane (Fig. 16). It is noted that the mesh is taken sufficiently large in order that the molten pool does not reach the boundary of the domain computation. This makes it possible to avoid the problems related to the definition of the normal at nodes of the mesh boundaries. Only half of the structure is modeled, and a single-pass hybrid weld is considered with a welding speed equal to $v_{welding} = 6 \text{ mm} \cdot \text{s}^{-1}$. The mesh (Fig. 16) includes 77,603 linear tetrahedral elements (P1 + /P1), 10,550 linear triangular elements (skin elements), and 14,922 nodes. In the fluid zone, the element size is 500 µm, and it can reach 5000 µm in the solid zone. The multiphasic ferritic steel grade S355J2G3 is used. All the data pertaining to this material are given in [25].

To study the effect of the fluid flow on the temperature field and the weld pool morphology, two simulations are carried out. The first simulation is a thermal-fluid computation (with fluid flow). Here, the heating time is 1.4 s, including 1 s for the creation of the molten pool (thermal computation with static heat sources and without fluid flow), and 0.4 s for the thermal-fluid simulation (with fluid flow) with the motion of heat sources. The time step is equal to 0.1 s for the creation of the molten pool and $2.5 \cdot 10^{-4} \text{ s}$ after the activation of the fluid flow. A temperature-dependent surface tension $\gamma = \gamma_L + \alpha(T - T_L)$ with $\gamma_L = 0.15 \text{ N} \cdot \text{m}^{-1}$, $T_L = 1500^{\circ}\text{C}$, and $\alpha = 10^{-4} \text{ N} \cdot \text{m}^{-1} \cdot \text{k}^{-1}$ is considered to include both the normal and Marangoni effects. The second simulation is a purely thermal simulation (without fluid flow) during 1.4 s.

In both simulations, heat exchanges due to convection and radiation on the boundary with the external medium are taken into account:

$$ \emptyset = \emptyset_{\text{convection}} + \emptyset_{\text{radiation}} \qquad (34) $$

with

$$ \emptyset_{\text{convection}} = h(T - T_0) \qquad (35) $$

$$ \emptyset_{\text{radiation}} = \epsilon \sigma_0 (T^4 - T_0^4) \qquad (36) $$

where $h = 40 \text{ W} \cdot \text{m}^{-2} \cdot \text{K}^{-1}$ is the transfer coefficient, $T_0 = 293 \text{ K}$ the room temperature, $\sigma_0$ the Stefan–Boltzmann constant, and $\epsilon = 0.7$ designates the emissivity. These conditions are enforced on all mesh faces except for the symmetry plane (YZ), where an insulating condition is imposed.

## 4.2. Results and discussion

The thermal-fluid simulation gives the distribution of the thermal field and fluid flow during the welding processes. In this simulation, the solid is schematized as a fluid with a high dynamic viscosity $\mu = 10^6 \text{ kg} \cdot \text{m}^{-1} \cdot \text{s}^{-1}$. The fluid is supposed to be Newtonian, with a dynamic viscosity $\mu = 2.5 \cdot 10^{-3} \text{ kg} \cdot \text{m}^{-1} \cdot \text{s}^{-1}$. The fluid flow is taken into account through the surface tension (normal and tangential effects) and the buoyancy. The fluid flow is activated at $t = 1 \text{ s}$, the time necessary to form a sufficient quantity of fluid in the weld pool.

In this section, a comparison of the results of the purely thermal and the thermal-fluid simulations are given (Fig. 17). Then, some illustrations of the results of the thermal-fluid simulation are given (Figs. 18–21). Fig. 17 illustrates the field
<figure><img src="image_20.png" /><figcaption>Fig. 17. Temperature field at $t = 1.12$ s (YZ plane).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 18. Temperature field, and velocity vector of the fluid flow at $t = 1.02$ s (YZ plane).</figcaption></figure>

temperature for purely thermal and thermal-fluid computations at $t = 1.12$ s. Figs. 18 and 19 show respectively the temperature field and the velocity field of the fluid flow at two instants ($t = 1.02$ s, $t = 1.13$ s), the first when the heat sources are fixed, and the second during the motion of the heat sources. Figs. 20 and 21 illustrate the evolution of the free surface in the weld pool at the same times ($t = 1.02$ s, $t = 1.13$ s). Here, the deformation of the free surface is induced by the surface tension and the buoyancy in the molten pool.

The effect of the fluid flow on the temperature field is very significant (Fig. 17). With regard to the thermal-fluid results, we note a significant temperature drop. This is explained by the positive derivative of the surface tension with respect to temperature which generates a penetrating weld pool, bringing the cold material from the edge towards the center of the weld pool.

During the thermal-fluid simulation, an average velocity approximately equal to $0.3 \ \text{m} \cdot \text{s}^{-1}$ is found. This value is consistent with literature values [56,57], whereas a maximal value of velocity equal to $0.6 \ \text{m} \cdot \text{s}^{-1}$ is noted when the heat sources begin to move. This increase of velocity is linked to the high-temperature gradient. Indeed, the gradient of the surface tension is proportional to the temperature gradient: therefore, when the temperature gradient increases, the velocity increases. Also, the higher $\frac{\partial \gamma}{\partial T}$, the higher the Marangoni number and the velocity of the fluid flow. The two effects of the surface tension with a positive $\frac{\partial \gamma}{\partial T}$ are clearly visible in Figs. 18 and 19. These figures show that the fluid flow moves from the edge of the molten zone to its center. Thus, a penetrating vortex is produced. This vortex is responsible for the change in the weld pool morphology. Tables 2 and 3 show the weld pool morphology at different times of the “purely” thermal and thermal-fluid computations. Maximal gaps equal to 15% for the length, 9.6% for the height and 1.6% for the width of the weld pool are found. A positive $\frac{\partial \gamma}{\partial T}$ gives a deeper and narrower weld pool. This effect of the fluid flow on the temperature
<figure><img src="image_22.png" /><figcaption>Fig. 19. Temperature field, and velocity vector of the fluid flow at $t = 1.13$ s (YZ plane).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 20. Vertical displacement (Uz) of the free surface at $t = 1.02$ s (YZ plane).</figcaption></figure>

field and the weld pool morphology can influence the microstructure and the mechanical computations. During a simulation taking into account only the buoyancy effect (without the tension surface effects), a maximal value equal to $0.03 \text{ m} \cdot \text{s}^{-1}$ is found. Therefore, it is about ten times smaller than that generated by the surface tension. This shows that the fluid flow linked to the surface tension specially the tangential effect "Marangoni effect" is dominant in this kind of process. Regarding the free surface, the maximal value of the vertical displacement is 0.5 mm, when the heat sources are fixed (Fig. 20), and it reaches 0.8 mm during the motion of the sources (Fig. 21). These values are similar to the literature results [18,58].
<figure><img src="image_24.png" /><figcaption>Fig. 21. Vertical displacement (Uz) of the free surface at $t = 1.13$ s (YZ plane).</figcaption></figure>

<table><caption>**Table 2**

Weld pool morphology at $t = 1.05$ s.</caption><thead><tr><th>[mm]</th><th>Thermal-fluid simulation</th><th>Purely thermal simulation</th><th>Gap</th></tr></thead><tbody><tr><td>Length</td><td>10.85</td><td>10.88</td><td>0.3%</td></tr><tr><td>Width</td><td>4.78</td><td>4.86</td><td>1.6%</td></tr><tr><td>Height</td><td>2</td><td>1.93</td><td>3.5%</td></tr></tbody></table>

<table><caption>**Table 3**

Weld pool morphology at $t = 1.4$ s.</caption><thead><tr><th>[mm]</th><th>Thermal-fluid simulation</th><th>Purely thermal simulation</th><th>Gap</th></tr></thead><tbody><tr><td>Length</td><td>13.7</td><td>16.14</td><td>15%</td></tr><tr><td>Width</td><td>4.78</td><td>4.86</td><td>1.6%</td></tr><tr><td>Height</td><td>3</td><td>2.71</td><td>9.6%</td></tr></tbody></table>

## 5. Conclusion

In order to study the interaction between the fluid flow and the solid deformation during the numerical simulation of the welding process, a new transient formulation for the numerical modeling of a weld pool was developed. It consists in considering the fluid flow in the weld pool through the surface tension and the buoyancy. In this context, a new method was used to simulate the surface tension. It incorporates both the normal force exerted onto the free surface ("curvature effect") and the tangential effect ("Marangoni effect") in a very simple and efficient way compared with classical methods. However, taking into account the normal effect of the surface tension still requires the simulation of the free surface. That is the reason why a new ALE approach was developed for modeling the evolution of the free surface. This method allows taking into account the solid deformation in a simple way, which is not the case for Eulerian methods.

This work was validated using two numerical examples. The first one is the sloshing test; it consists in simulating the oscillations of a liquid in a container, and it allowed validating the normal effect of the surface tension and the ALE approach. The second one is the benchmark problem, it represents the numerical simulation of spot welding "GTWA", and it allowed validating the tangential effect of the surface tension.

After validation, the 3D thermal-fluid simulation of a hybrid laser/arc was carried out. It gives good results in terms of the velocity of the fluid flow and the evolution of the free surface compared to the literature results. To understand the effect of the fluid flow during the simulation of hybrid welding, a comparison of results between purely thermal and thermal-fluid simulations was performed. It shows the significant influence of the fluid flow on the temperature field and the morphology of the weld pool. For example, we note a significant temperature drop in the thermal-fluid simulation. This may be explained by the positive derivative of the surface tension with respect to temperature, which brings the cold
material from the edge towards the center of the weld pool. Regarding the molten pool morphology, maximal gaps equal to 15% for the length, 9.6% for the height, and 1.6% for the width of the weld pool were found. The fluid flow induced by a positive derivative of the surface tension gives a deeper and narrower weld pool. The effect of the fluid flow on the temperature field and the weld pool morphology can influence the microstructure and the mechanical fields of the solid part and, on the counterpart, the deformations of the solid part can modify the fluid flow in the molten pool. A work is in progress to couple, in an ALE formulation, the fluid flow in the molten pool with the deformations of the solid part. The solid/fluid interaction will be implemented by imposing solid velocities during the fluid computation, the solid velocities coming at each time step from a solid computation.

## References

[1] J.-M. Bergheau, J.-B. Leblond, Coupling between heat flow, metallurgy and stress-strain computations in steels – the approach developed in the computer code SYSWELD for welding or quenching, in: M. Rappaz, M.R. Ozgu, K.W. Mahin (Eds.), Modeling of Casting, Welding and Advanced Solidification Processes V, The Minerals, Metals & Materials Society, 1991, pp. 203–210.

[2] P. Duranton, J. Devaux, V. Robin, P. Gilles, J.-M. Bergheau, 3D modelling of multipass welding of a 316L stainless steel pipe, J. Mater. Process. Technol. 153-154 (2004) 457–463.

[3] M. Zain ul Abdein, D. Nelias, J.F. Julien, D. Deloison, Prediction of laser beam welding-induced distortions and residual stresses by numerical simulation for aeronautic application, J. Mater. Process. Technol. 209 (6) (2009) 2907–2917.

[4] P. Martinson, S. Daneshpour, M. Kocak, S. Riekehr, P. Staron, Residual stress analysis of laser spot welding of steel sheets, Mater. Des. 30 (9) (2009) 3351–3359.

[5] É. Feulvarch, V. Robin, J.-M. Bergheau, Thermometallurgical and mechanical modelling of welding – application to multipass dissimilar metal girth welds, Sci. Technol. Weld. Join. 16 (3) (2011) 221–226.

[6] S.A. Tsirkas, P. Papanikos, K. Pericleous, N. Strusevich, F. Boitout, J.-M. Bergheau, Evaluation of distortions in laser welded shipbuilding parts using a local-global finite element approach, Sci. Technol. Weld. Join. 8 (2) (2003) 79–88.

[7] M. Zain ul Abdein, D. Nelias, J.F. Julien, D. Deloison, Experimental investigation and finite element simulation of laser beam welding induced residual stresses and distortions in thin sheets of (AA) 6056-t4, Mater. Sci. Eng. A 527 (12) (2010) 3025–3039.

[8] Z. Zhang, P. Ge, G.Z. Zhao, Numerical studies of post weld heat treatment on residual stresses in welded impeller, Int. J. Press. Vessels Piping 153 (2017) 1–14.

[9] J.-D. Caprace, G. Fu, J.F. Carra, H. Remes, S.B. Shin, A benchmark study of uncertainty in welding simulation, Mar. Struct. 56 (2017) 69–84.

[10] G.M.L. Xiong, C. Wang, X. Hu, Y. Wei, A thermal-metallurgical-mechanical model for laser welding Q235 steel, J. Mater. Process. Technol. 238 (2016) 39–48.

[11] J. Cheon, S.-J. Na, Prediction of welding residual stress with real-time phase transformation by CFD thermal analysis, Int. J. Mech. Sci. 131–132 (2017) 37–51.

[12] H. Zhao, W. Niu, B. Zhang, Y. Lei, M. Kodama, T. Ishide, Modelling of keyhole dynamics and porosity formation considering the adaptive keyhole shape and three-phase coupling during deep-penetration laser welding, J. Phys. D, Appl. Phys. (2011), https://doi.org/10.1088/0022-3727/44/48/485302.

[13] C. Panwisawas, Y. Sovani, R.P. Turner, J.W. Brooks, H.C. Basoalto, I. Choquet, Modelling of thermal fluid dynamics for fusion welding, J. Mater. Process. Technol. (2017), https://doi.org/10.1016/j.jmatprotec.2017.09.019.

[14] C-c. Liu, J-s. He, Numerical analysis of thermal fluid transport behavior during electron beam welding of 2219 aluminum alloy plate, Trans. Nonferr. Met. Soc. China 27 (2017) 1319–1326.

[15] M.R. Nekouie Esfahani, J. Coupland, S. Marimuthu, Numerical simulation of alloy composition in dissimilar laser welding, J. Mater. Process. Technol. 224 (2015) 135–142.

[16] M. Bachmann, V. Avilov, A. Gumenyuk, M. Rethmeier, Numerical assessment and experimental verification of the influence of the Hartmann effect in laser beam welding processes by steady magnetic fields, Int. J. Therm. Sci. 101 (2016) 24–34.

[17] X. Meng, G. Qin, Z. Zou, Investigation of humping defect in high speed gas tungsten arc welding by numerical modelling, Mater. Des. 94 (2016) 69–78.

[18] M.C. Nguyen, M. Medale, O. Asserin, S. Gounand, P. Gilles, Sensitivity to welding positions and parameters in GTA welding with a 3D multiphysics numerical model, Numer. Heat Transf., Part A, Appl. 71 (3) (2017) 233–249.

[19] T.D. Anderson, J.N. DuPont, T. DebRoy, Origin of stray grain formation in single-crystal superalloy weld pools from heat transfer and fluid flow modeling, Acta Mater. 58 (2010) 1441–1454.

[20] H.L. Wei, J.W. Elmer, T. DebRoy, Three-dimensional modeling of grain structure evolution during welding of an aluminum alloy, Acta Mater. 126 (2017) 413–425.

[21] H.L. Wei, J.W. Elmer, T. DebRoy, Crystal growth during keyhole mode laser welding, Acta Mater. 126 (2017) 413–425.

[22] L.L.H.L. Wei, J.X. Xue, T. DebRoy, Special features of double pulsed gas metal arc welding, J. Mater. Process. Technol. 251 (2018) 369–375.

[23] H. Sallem, T. Sayet, É. Feulvarch, J.-B. Leblond, J.-M. Bergheau, A comparison of two modeling approaches for the prediction of residual stresses caused by the welding process, in: NUMIFORM, 2016.

[24] Z. Zhuanzhuan, W. ChuanSong, Effect of fluid flow in the weld pool on the numerical simulation accuracy of the thermal field in hybrid welding, J. Manuf. Process. 20 (2015) 215–223.

[25] H. Amin-El-Sayed, Simulation numérique du soudage : couplage des écoulements dans le bain fondu avec les déformations de la partie solide, PhD thesis, University of Lyon, France, 2014.

[26] Y. Saadlaoui, H. Sallem, É. Feulvarch, T. Sayet, J.-B. Leblond, J.-M. Bergheau, Influence des écoulements de la matière sur l'évolution des contraintes résiduelles durant le procédé de soudage, in: 13ᵉ colloque national en calcul des structures CSMA, 15–19 May 2017, Giens, France, 2017.

[27] K.J. Ruschak, A method for incorporating free surface boundaries with surface tension in finite element fluid-flow simulators, Int. J. Numer. Methods Fluids 15 (1980) 639–648.

[28] P. Bach, O. Hassager, An algorithm for the use of the Lagrangian specification in Newtonian fluid mechanics with applications to free surface fluid flow, J. Fluid Mech. 152 (1985) 173–190.

[29] K.N. Christodoulou, L.E. Scriven, The fluid mechanics of slide coating, J. Fluid Mech. 208 (1989) 321–354.

[30] R.A. Cairncross, P.R. Schunk, T.A. Baer, R.R. Bao, P.A. Sackinger, Finite element method for free surface flows of incompressible fluids in three dimensions. Part I. Boundary fitted mesh motion, Int. J. Numer. Methods Fluids 33 (2000) 375–403.

[31] M. Bellet, Implementation of surface tension with wall adhesion effects in a three-dimensional finite element model for fluid flow, Commun. Numer. Methods Eng. 17 (2001) 563–579.

[32] S. Rabier, M. Medale, Computation of free surface flows with a projection FEM in a moving mesh framework, Comput. Methods Appl. Mech. Eng. 192 (2003) 4703–4721.

[33] C.W. Hirt, A.A. Amsden, J.L. Cook, An arbitrary Lagrangian–Eulerian computing method for all flow speeds, J. Comput. Phys. 135 (1997) 203–216.
[34] C.W. Hirt, B.D. Nichols, Volume of fluid (vof) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–225.

[35] S. Osher, J.A. Sethian, Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49.

[36] F. Lu, X. Li, Z. Li, X. Tang, H. Cui, Formation and influence mechanism of keyhole-induced porosity in deep-penetration laser welding based on 3D transient modeling, Int. J. Heat Mass Transf. 90 (2015) 1143–1152.

[37] S. Pang, W. Chen, W. Wang, A quantitative model of keyhole in stability induced porosity in laser welding of titanium alloy, Metall. Mater. Trans. A 45 (6) (2014) 2808–2818.

[38] M. Medale, S. Rabier, C. Xhaard, A thermo-hydraulic numerical model for high energy welding processes, Rev. Eur. Éléments Finis 13 (3-4) (2004) 207–229.

[39] M. Medale, C. Touvrey, R. Fabbro, An axi-symmetric thermo-hydraulic model to better understand spot laser welding, Eur. J. Mech. 17 (5-7) (2008) 795–806.

[40] M. Medale, C. Xhaard, R. Fabbro, A thermo-hydraulic numerical model to study spot laser welding, C. R. Mecanique 335 (2007) 280–286.

[41] J.-B. Leblond, H. Amin El Sayed, J.-M. Bergheau, On the incorporation of surface tension in finite-element calculations, C. R. Mecanique 341 (2013) 770–775.

[42] J.-M. Bergheau, R. Fortunier, Finite Element Simulation of Heat Transfer, ISTE-Wiley, ISBN 978-1-84821-053-0, 2008, 279 pages.

[43] É. Feulvarch, J.-M. Bergheau, J.-B. Leblond, An implicit finite element algorithm for the simulation of diffusion with phase changes in solids, Int. J. Numer. Methods Eng. 78 (2009) 1492–1512.

[44] W. Dettmer, D. Peric, A computational framework for free surface fluid flows accounting for surface tension, Comput. Methods Appl. Mech. Eng. 195 (2006) 3038–3071.

[45] P.H. Saksono, D. Peric, On finite element modelling of surface tension. Variational formulation and applications – Part I: quasistatic problems, Comput. Mech. 38 (2006) 265–281.

[46] Y. Demay, A. Béliveau, A. Fortin, A numerical method for the deformation of two-dimensional drops with surface tension, Int. J. Comput. Fluid Dyn. 10 (1998) 225–240.

[47] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, G. Zanetti, Modelling merging and fragmentation in multiphase flows with SURFER, J. Comput. Phys. 113 (1) (1994) 134–147.

[48] T. Heuzé, J.-B. Leblond, J.-M. Bergheau, Modélisation des couplages fluide/solide dans les procédés d'assemblage à haute température, Méc. Ind. 12 (2011) 183–191.

[49] B. Ramaswamy, Numerical simulation of unsteady viscous free surface flow, J. Comput. Phys. 90 (1990) 396–430.

[50] S. Popinet, S. Zaleski, A front-tracking algorithm for accurate representation of surface tension, Int. J. Numer. Methods Fluids 30 (1999) 775–793.

[51] A. Prosperetti, Motion of two superposed viscous fluids, Phys. Fluids 24 (7) (1981) 1217–1223.

[52] P. Girard, M. Bellet, G. Caillibote, M. Carin, S. Gounant, F. Mathey, M. Médale, Benchmark for Fluid Flow in Weld Pool Simulation: Two Dimensional Transient Computational Models for Arc Welding, Internal research report, 2005.

[53] M. Hamide, Modélisation numérique du soudage à l'arc des aciers, PhD thesis, École des mines de Paris, 2008.

[54] ESI Group, User Manual, 2015.

[55] J. Goldack, A. Chakravati, M. Bibby, Metall. Trans. 15B (1984) 299–305.

[56] X.S. Gao, C.S. Wu, S.F. Goecke, H. Kugler, Numerical simulation of temperature field, fluid flow and weld bead formation in oscillating single mode laser-GMA hybrid welding, J. Mater. Process. Technol. 242 (2017) 147–159.

[57] B. Huang, X. Chen, S. Pang, R. Hu, A three-dimensional model of coupling dynamics of Keyhole and weld pool during electron beam welding, Int. J. Heat Mass Transf. 115 (2017) 159–173.

[58] M.C. Nguyen, M. Medale, O. Asserin, S. Gounand, Modélisation et simulation multiphysique du bain de fusion en soudage à l'arc TIG, PhD thesis, University of Aix-Marseille, France, 2015.