Numerical study of molten metal melt pool behaviour during conduction-mode laser spot melting

Laser spot melting of a metallic S705 alloy, as shown schematically in figure 1 and as experimentally studied by Pitscheneder et al [18], was numerically simulated as a representative example in the present work. A defocused laser-beam with a radius of $r_\mathrm{b} = 1.4 \, \textrm{mm}$ heats the bulk material from its top surface for $4 \, \textrm{s}$. The laser-beam power Q is set to $3850 \, \textrm{W}$ with a top-hat intensity distribution. The averaged laser absorptivity of the material surface η is assumed to remain constant at 13% [18]. The absorbed laser power leads to an increase in temperature and subsequent melting of the base material. The base material is a rectangular cuboid shape with a base size (L × L) of $24 \times 24 \, \textrm{mm}^2$ and a height ($H_\mathrm{m}$) of $10 \, \textrm{mm}$, initially at an ambient temperature of $T_\mathrm{i} = 300~\textrm{K}$. A layer of air with a thickness of $H_\mathrm{a} = 2~\textrm{mm}$ is considered above the base material to monitor the gas-metal interface evolution. Except for the surface tension, the material properties are assumed to be constant and temperature independent and are presented in table 1. Although temperature-dependent properties can be employed in the present model to enhance the model accuracy, employing temperature-independent properties facilitates comparison of the results of the present work with previous results published in the literature [10, 19, 20], where temperature-independent material properties are employed. Moreover, further studies are required to enhance the accuracy of calculation and measurement of temperature-dependent material properties, particularly for the liquid phase above the melting temperature. The effects of employing temperature-dependent material properties on the thermal and fluid flow fields and melt-pool shape are discussed in detail in section 4.3. A change in surface-tension due to the non-uniform temperature distribution over the gas–liquid interface induces thermocapillary stresses that drive the melt flow. This fluid motion from low to high surface-tension regions changes the temperature distribution in the melt pool [21] and can lead to surface deformations that change the melt-pool shape and properties of the material after solidification [22]. The sulphur contained in the alloy can alter the surface tension of the molten material and its variations with temperature [23].

Zoom In Zoom Out Reset image size

Table 1. Thermophysical properties of the Fe-S alloy and air used in the present study. Values are taken from [20].

Property	Fe-S alloy	Air	Unit
Density ρ	8100	1.225	kg m−3
Specific heat capacity $c_\mathrm{p}$	670	1006	J kg−1 K−1
Thermal conductivity k	22.9	0.024	W m−1 K
Viscosity µ	6 × 10−3	1.8 × 10−05	kg m−1 s−1
Latent heat of fusion $L_\mathrm{f}$	250800	–	J kg−1
Liquidus temperature $T_\mathrm{l}$	1620	–	K
Solidus temperature $T_\mathrm{s}$	1610	–	K

A three-dimensional multiphase model based on the finite-volume method was utilised to predict the melt pool dynamic behaviour and the associated transport phenomena. The present model is developed for conduction-mode laser melting. Further considerations will be required to develop the model for keyhole mode welding, including the complex laser-matter interactions, changes in surface chemistry and consequent surface tension variations for non- or partially shielded welding conditions. Both the molten metal and air were treated as Newtonian and incompressible fluids. Thermal buoyancy forces, which in laser spot melting are generally negligible compared to thermocapillary forces [24], were neglected both in the liquid metal and in the air. During conduction-mode laser-melting, the heat-source energy density is too low to cause significant vaporisation of the liquid metal [25], and the surface deformations are too small to cause multiple reflections of the laser rays at the liquid surface [26]. Furthermore, our preliminary results indicate that the maximum temperature of the melt pool is predominantly below $2700 \, \textrm{K}$, significantly less than the boiling temperature of stainless steels, which is typically $\mathcal{O}(3100) \, \mathrm{K}$. Hence, multiple reflections and vaporisation were ignored. With this, the unsteady conservation equations of mass, momentum and energy were defined as follows:

$$\begin{equation} \nabla\cdot\left(\rho\vec{V}\right) = 0, \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial}{\partial t} \left(\rho \vec{V}\right) + \nabla\cdot\left(\rho\vec{V} \vec{V}\right) = -\nabla p + \nabla\cdot\left(\mu\nabla\vec{V}\right) + \vec{S}_\mathrm{d} + \vec{F}_\mathrm{s}, \end{equation} \tag{ 2 }$$

$$\begin{align} \frac{\partial}{\partial t} \left(\rho h\right) + \nabla\cdot\left(\rho\vec{V} h\right) & = \nabla\cdot\left(\frac{k}{c_\mathrm{p}}\nabla h\right) - \frac{\partial}{\partial t} \left(\rho \Delta H\right)\nonumber\\ &\quad - \nabla\cdot\left(\rho\vec{V}\Delta H\right) + S_\mathrm{T}, \end{align} \tag{ 3 }$$

where ρ is density, t time, $\vec{V}$ velocity vector, p pressure, µ dynamic viscosity, h sensible heat, k thermal conductivity, $c_\mathrm{p}$ specific heat capacity at constant pressure and ΔH latent heat. The enthalpy of the material H was defined as the sum of the latent heat and the sensible heat [27], and is expressed as

$$\begin{equation} H = h + \Delta H = \left(h_\mathrm{ref} + \int_{T_\mathrm{ref}}^{T} c_\mathrm{p} \mathop{}\mathrm{d}T\right) + f_\mathrm{L}L_\mathrm{f}, \end{equation} \tag{ 4 }$$

where $h_\mathrm{ref}$ is reference enthalpy, $T_\mathrm{ref}$ reference temperature, $L_\mathrm{f}$ latent heat of fusion, and $f_\mathrm{L}$ local liquid volume-fraction. Assuming the liquid volume-fraction in the metal to be a function of temperature only, a linear function [27] was used to calculate the liquid volume-fraction as follows:

$$\begin{equation} f_\mathrm{L} = \frac{T - T_\mathrm{s}}{T_\mathrm{l} - T_\mathrm{s}}; \quad T_\mathrm{s} \le T \le T_\mathrm{l}, \end{equation} \tag{ 5 }$$

where $T_\mathrm{s}$ and $T_\mathrm{l}$ are the solidus and liquidus temperatures, respectively.

To suppress fluid velocities in solid regions and to model fluid flow damping in the so-called "mushy zone", where phase transformation occurs within a temperature range between $T_\mathrm{s}$ and $T_\mathrm{l}$, a momentum sink term $\vec{S}_\mathrm{d}$ [28] was implemented into the momentum equation as

$$\begin{equation} \vec{S}_\mathrm{d} = C\ \frac{(1 - f_\mathrm{L})^2}{f_\mathrm{L}^3 + \epsilon} \ \vec{V}, \end{equation} \tag{ 6 }$$

where C is the mushy-zone constant and is a small number to avoid division by zero as $f_\mathrm{L}$ approaches 0. The mushy-zone constant C was set to $10^7\,\textrm{kg}\, \textrm{m}^{-2}\, \textrm{s}^{-2}$, which is large enough to dampen fluid velocities in solid regions according to the criterion defined by Ebrahimi et al [29], and was equal to 10⁻³.

To capture the position of the gas–metal interface during melting, the volume-of-fluid (VOF) method developed by Hirt and Nichols [14] was utilised. This requires the solution of a transport equation for one additional scalar variable φ (i.e., the so-called volume-of-fluid fraction) that varies between 0 in the gas phase and 1 in the metal phase:

$$\begin{equation} \frac{\partial \phi}{\partial t} + \nabla\cdot\left(\phi\vec{V}\right) = 0. \end{equation} \tag{ 7 }$$

Computational cells with 0 ≤ φ ≤ 1 are in the gas-metal interface region. The effective thermophysical properties of the fluid is then calculated using the mixture model as follows [30]:

$$\begin{equation} \zeta = \phi \, \zeta_\mathrm{metal} + \left(1-\phi\right) \zeta_\mathrm{gas}, \end{equation} \tag{ 8 }$$

where ζ corresponds to density ρ, viscosity µ, thermal conductivity k and specific heat capacity $c_\mathrm{p}$.

Based on the method developed by Brackbill et  al [31], surface tension and thermocapillary forces acting on the gas-metal interface were applied as volumetric forces in these interface cells and were introduced into the momentum equation as a source term $\vec{F}_\mathrm{s}$ as follows:

$$\begin{equation} \vec{F}_\mathrm{s} = \vec{f}_\mathrm{s} \lVert \nabla\phi \rVert\frac{2\rho}{\rho_\mathrm{gas} + \rho_\mathrm{metal}}, \end{equation} \tag{ 9 }$$

where subscripts indicate the phase. The term $2\rho/\!\left(\rho_\mathrm{gas} + \rho_\mathrm{metal}\right)$ was utilised to abate the effect of the large metal-to-gas density ratio by redistributing the volumetric surface-forces towards the metal phase (i.e., the heavier phase). $\vec{f}_\mathrm{s}$ is the surface force per unit area and was defined as follows:

$$\begin{equation} \vec{f}_\mathrm{s} = \sigma\kappa\vec{n} + \frac{\mathrm{\textrm{d}}\sigma}{\mathrm{d}T} \left[\nabla T - \vec{n}\left(\vec{n}\cdot\nabla T\right)\right], \end{equation} \tag{ 10 }$$

where the first term on the right-hand side is the normal component, which depends on the value of surface tension and curvature of the melt-pool surface, and the second term is the tangential component that is affected by the surface tension gradient. In equation (10), σ is surface tension, and κ the surface curvature is

$$\begin{equation} \kappa = -\left(\nabla\cdot\vec{n}\right), \end{equation} \tag{ 11 }$$

and $\vec{n}$ the surface unit normal vector is

$$\begin{equation} \vec{n} = \frac{\nabla\phi}{\lVert \nabla \phi \rVert}. \end{equation} \tag{ 12 }$$

2.2.1. Surface tension.

The temperature dependence of the surface tension of a liquid solution with a low concentration of surfactant can be approximated using a theoretical correlation derived on the basis of the combination of Gibbs and Langmuir adsorption isotherms as follows [32, 33]:

$$\begin{equation} \sigma = \sigma^\circ - {R}\, T\, \Gamma_\mathrm{s}\, \ln(1 + K a_\mathrm{s}), \end{equation} \tag{ 13 }$$

where σ is the surface tension of the solution, σ° the pure solvent surface tension, R the gas constant, $\Gamma_\mathrm{s}$ the adsorption at saturation, K the adsorption coefficient, and $a_\mathrm{s}$ the activity of the solute. From this, Sahoo et al [23] derived a correlation for binary molten metal–surfactant systems, including binary Fe–S systems:

$$\begin{align} \sigma & = \sigma_{m}^\circ + \left(\frac{\partial\sigma}{\partial T}\right)^\circ \left(T - T_{m}\right)\nonumber\\ & \quad - {R}\, T\, \Gamma_\mathrm{s}\, \ln\left[1 + \psi\, a_\mathrm{s} \exp\left(\frac{-\Delta H^\circ}{\mathrm{R} T}\right)\right], \end{align} \tag{ 14 }$$

where $\sigma_{m}^\circ$ is the surface tension of pure molten-metal at the melting temperature T_m , ψ an entropy factor, ΔH° the standard heat of adsorption, and $\left(\partial \sigma / \partial T\right)^\circ$ the temperature coefficient of the surface tension of the pure molten-metal. Values of the properties used in equation (14) to calculate the temperature dependence of the surface tension of the molten Fe–S alloy are presented in table 2. Variations of the surface tension and its temperature coefficient are shown in figure 2 for an Fe–S alloy with 150 ppm sulphur content.

Zoom In Zoom Out Reset image size

Table 2. The values used to calculate the surface tension of molten Fe–S alloy [23]

Property	Value	Unit
$\sigma_\mathrm{m}^\circ$	1.943	N m−1
$\left(\partial \sigma / \partial T\right)^\circ$	−4.3 × 10−4	N m−1 K−1
$\Gamma_\mathrm{s}$	1.3 × 10−5	mol m−2
ψ	3.18 × 10−3	–
$a_\mathrm{s}$	150	–
ΔH°	−166.2 × 105	J mol−1

2.2.2. Laser heat source.

The laser supplies a certain amount of energy (η  Q) to the material. Its energy input was modelled by adding a volumetric source term $S_\mathrm{T}$ to the energy equation (equation (3)) in cells at the gas–metal interface. The surface deformations affect the total energy input to the material [34]. This effect is generally not considered in numerical simulations of melt pool behaviour. To investigate the influence of surface deformations, five different methods for the heat-source implementation were considered.

2.2.2.1. Case 1, heat-source implementation without taking the influence of surface deformations on total energy input into consideration.

In this case, the heat source model is defined as

$$\begin{equation} S_\mathrm{T} = \left\{\!\!\!\! \begin{array}{ll} \frac{\eta\, Q}{\pi r_\mathrm{b}^2} \lVert \nabla\phi \rVert\frac{2\,\rho\, c_\mathrm{p}}{\left(\rho\, c_\mathrm{p}\right)_\mathrm{gas} + \left(\rho\, c_\mathrm{p}\right)_\mathrm{metal}}, & \textrm{if}\,r \le r_\mathrm{b} \\ 0, & \text{otherwise} \end{array}\right. \end{equation} \tag{ 15 }$$

where r is radius defined as $\sqrt{x^2 + y^2}$. This method is the most common method in modelling melt-pool surface oscillations without having a deep penetration (keyhole formation) (see for instance, [25, 26]). However, variations in the total energy input caused by a change in melt-pool surface shape, through changes in $\lVert \nabla\phi \rVert$, causing $\iiint \limits_{\mathrm{\forall}} S_\mathrm{T} \mathop{}\!\mathrm{d}V$ to be different from ηQ, are an inherent consequence of using this method. Additionally, every segment of the melt-pool surface is exposed to the same amount of energy input regardless of the local surface orientation.

2.2.2.2. Case 2, heat source adjustment.

To conserve the total energy input, an adjustment coefficient ξ was introduced to the heat source model as follows:

$$\begin{equation} S_\mathrm{T} = \left\{\!\!\!\! \begin{array}{ll} \frac{\xi\, \eta\, Q}{\pi r_\mathrm{b}^2} \lVert \nabla\phi \rVert\frac{2\,\rho\, c_\mathrm{p}}{\left(\rho\, c_\mathrm{p}\right)_\mathrm{gas} + \left(\rho\, c_\mathrm{p}\right)_\mathrm{metal}}, & \textrm{if}\,r \le r_\mathrm{b} \\ 0, & \text{otherwise} \end{array}\right. \end{equation} \tag{ 16 }$$

where ξ was evaluated at the beginning of every time-step and was defined as

$$\begin{equation} \xi = \frac{\eta \, Q}{\iiint \limits_{\mathrm{\forall}} S_\mathrm{T} \mathop{}\mathrm{d}V}, \end{equation} \tag{ 17 }$$

where ∀ stands for the computational domain. Utilising this technique for welding simulations has already been reported in the literature (see for instance, [15, 16, 34, 35]). Variations in energy absorption with deformations of the melt-pool surface are, however, not taken into consideration in this method.

2.2.2.3. Case 3, heat source redistribution.

In this case, the energy flux was redistributed over the melt-pool surface assuming the energy absorption is a function of the local surface orientation. Accordingly, the local energy absorption is a maximum where the surface is perpendicular to the laser ray and is a minimum where the surface is aligned with the laser ray [36, 37]. This represents a simplified absorptivity model based on the Fresnel's equation [38] and is expressed mathematically, assuming the laser rays being parallel and in the z-direction, as follows:

$$\begin{align} S_\mathrm{T} = \left\{ \begin{array}{ll} \frac{\eta\, Q}{\pi r_\mathrm{b}^2} \lVert \nabla\phi \rVert\frac{2\,\rho\, c_\mathrm{p}}{\left(\rho\, c_\mathrm{p}\right)_\mathrm{gas} + \left(\rho\, c_\mathrm{p}\right)_\mathrm{metal}} \left| \frac{\nabla \phi \left[z\right]}{\lVert \nabla \phi \rVert} \right|. & \textrm{if}\, r \le r_\mathrm{b} \\ 0. & \text{otherwise} \end{array}\right. \end{align} \tag{ 18 }$$

the total energy input is not necessarily conserved in this method.

2.2.2.4. Case 4, heat source redistribution and adjustment.

Utilising the same technique introduced in Case 2, the heat source model in Case 3 was adjusted to guarantee that total energy input is conserved. Hence, the heat source model was defined as

$$\begin{equation} S_\mathrm{T} = \left\{\!\!\!\! \begin{array}{ll} \frac{\xi\, \eta\, Q}{\pi r_\mathrm{b}^2} \lVert \nabla\phi \rVert\frac{2\,\rho\, c_\mathrm{p}}{\left(\rho\, c_\mathrm{p}\right)_\mathrm{gas} + \left(\rho\, c_\mathrm{p}\right)_\mathrm{metal}} \left| \frac{\nabla \phi \left[z\right]}{\lVert \nabla \phi \rVert} \right|. & \textrm{if}\,r \le r_\mathrm{b} \\ 0. & \text{otherwise} \end{array}\right. \end{equation} \tag{ 19 }$$

2.2.2.5. Case 5, Flat non-deformable free surface.

In this case, the surface was assumed to remain flat. Thermocapillary shear stresses caused by a non-uniform temperature distribution on the gas-metal interface were applied as a boundary condition, hence modelling the gas phase was not required. Consequently, the total energy absorbed by the surface was fixed in this method. The boundary conditions are described in section 2.2.3.

Table 3 presents a summary of the cases considered in the present work.

Table 3. Summary of the cases studied in the present work and the features included in the model for each case.

Name	Deformable free-surface	Heat source adjustment	Heat source redistribution
Case 1	Yes	No	No
Case 2	Yes	Yes	No
Case 3	Yes	No	Yes
Case 4	Yes	Yes	Yes
Case 5	No	No	No

2.2.3. Boundary conditions.

The bottom and lateral surfaces of the metal part were modelled to be no slip walls, but in fact they remain solid during the simulation time. At the boundaries of the gas layer above the metal part, a constant atmospheric pressure was applied (i.e., $p = 101.325 \, \textrm{kPa}$) allowing air to flow in and out of the domain. The outer boundaries of the computational domain were modelled as adiabatic since heat losses through these boundaries are negligible compared to the laser power [20].

For Case 5, in which no gas layer is modelled explicitly, a thermocapillary shear-stress was applied as a boundary condition in the molten regions of the top surface. In the irradiated region on the material top-surface, a constant uniform heat flux was applied, while outside of this region the surface was assumed to be adiabatic [18]. The thermal and thermocapillary shear-stress boundary conditions were defined, respectively, as

$$\begin{equation} k \frac{\partial T}{\partial \vec{n}} = \left\{\!\!\!\! \begin{array}{ll} \frac{\eta\, Q}{\pi r_\mathrm{b}^2}, & \textrm{if}\,r \le r_\mathrm{b} \\ 0, & \text{otherwise} \end{array}\right. \end{equation} \tag{ 20 }$$

and

$$\begin{equation} -\mu \frac{\partial \vec{V_\mathrm{t}}}{\partial \vec{n}} = \frac{\mathop{}\mathrm{d}\sigma}{\mathop{}\mathrm{d}T} \frac{\partial T}{\partial \vec{\tau}}, \end{equation} \tag{ 21 }$$

where $\vec{V_\mathrm{t}}$ is the tangential velocity vector, and $\vec{\tau}$ the tangential vector to the top surface.

To verify the reliability and accuracy of the present numerical simulations, the melt-pool shapes obtained from the present simulations, for the problem introduced in section 2.1, are compared to experimental observations reported by Pitscheneder et  al [18]. Figure 3 shows a comparison between the melt-pool shape obtained from the numerical simulation after $5 \, \textrm{s}$ of heating and the post-solidification experimental observation, which indicates a reasonable agreement. The maximum absolute deviations between the present numerical predictions and experimental data for the melt-pool width and depth is less than 5% and 2%, respectively. However, it should be noted that the thermal conductivity and viscosity of the liquid material were artificially increased in the simulations by a factor $\mathcal{F} = 7$ with respect to their reported experimental values, as suggested by Pitscheneder et al [18]. This is known as employing an 'enhancement factor', mainly to achieve agreement between numerical and experimental results [19]. Independent studies conducted by Saldi et al [19] and Ehlen et al [43] revealed that without using such an enhancement factor, the fluid flow structure and the melt-pool shape in simulations differ drastically from experimental observations. The use of an enhancement factor is often justified by the possible occurrence of turbulence and its influence on heat and momentum transfer in the melt pool, which is assumed to be uniform in the melt pool. However, the high-fidelity numerical simulations conducted by Kidess et al [13, 20] on a melt pool with a flat non-deformable melt-pool surface revealed that turbulent enhancement is strongly non-uniform in the melt pool, resulting in an ω-shaped melt pool that differs notably from the results of simulations assuming uniform transport enhancement. We extended this study and performed a high-fidelity simulation based on the large eddy simulation (LES) turbulence model and took the effects of surface deformations into consideration [44] and showed that for the conduction-mode laser melting problem we considered, the influence of surface oscillations on melt pool behaviour is larger compared to the effects of turbulent flow in the melt pool; nevertheless the results do not agree with experiments without using some enhancement. These observations along with previous studies [19, 43, 45] suggest that the published weld pool models lack the inclusion of significant physics. The neglect in simulations of relevant physics such as chemical reactions and unsteady oxygen absorption by the melt-pool surface, non-uniform unsteady surfactant distribution over the melt-pool surface, re-solidification, free surface evolution and three-dimensionality of the fluid flow field has been postulated as reasons why such an ad hoc and unphysical enhancement factor is needed to obtain agreement with experimental data [46–48]. However, the inclusion of such factors will increase the model complexity and computational costs. Further investigations are essential to realising the problem with sufficient details. In section 4.4, the effects of the value used for the enhancement factor $\mathcal{F}$ on the thermal and fluid flow fields and melt-pool oscillatory behaviour are discussed in more detail.

Zoom In Zoom Out Reset image size

The experimental and numerical data reported by He et  al [49], who investigated unsteady heat and fluid flow in the melt pool during conduction-mode laser spot welding of stainless steel (AISI304) plates, were considered to validate the reliability of the present model. The melt pool shape obtained from the present model assuming a non-deformable free surface (Case 5) was compared with experimental and numerical data reported by He et al [49] and the results are shown in figure 4. In this problem, the laser power was set to $1967 \, \textrm{W}$, the beam radius was $570 \, \unicode{x03BC}\textrm{m}$ and laser pulse duration was $3 \, \textrm{ms}$. The thermal conductivity and viscosity of the molten material were artificially increased in the simulations by a factor $\mathcal{F} = 17$ and $\mathcal{F} = 11$, respectively, as suggested by He et al [49]. The maximum absolute deviations between the present numerical predictions and experimental data for the melt-pool width and depth is less than 8% and 2.5%, respectively. The results obtained from the present model are indeed in reasonable agreement with the reference data, indicating the validity of the present model in reproducing the results reported in literature.

Zoom In Zoom Out Reset image size

The reliability of the present model was also investigated by comparing the numerically predicted melt-pool size and surface deformations with experimental observations of Cunningham et al [50]. The evolutions of the melt-pool shape and its surface depression under stationary laser melting of a Ti-6Al-4V plate under conduction mode were studied. In this problem, the laser power is $156 \, \textrm{W}$ and has a Gaussian distribution and the radius of the laser spot is $7 \times 10^{-5} \, \textrm{m}$. The thermophysical properties of Ti-6Al-4V suggested by Sharma et  al [51] were employed for calculations. Since the melt pool surface temperature reaches the boiling point, the effects of recoil pressure $p_\mathrm{r}$ and evaporation heat loss $q_\mathrm{e}$ were included in the model. The recoil pressure [52] and evaporation heat loss [53] were determined respectively as follows:

$$\begin{equation} p_\mathrm{r} = 0.54 \cdot P_0 \exp\left(\frac{L_\mathrm{v} \, {M} \, \left(T - T_\mathrm{v}\right)}{{R} \, T \, T_\mathrm{v}}\right), \end{equation} \tag{ 22 }$$

$$\begin{equation} q_\mathrm{e} = -0.82 \frac{L_\mathrm{v} \, {M}}{\sqrt{2 \pi {M} \, {R} \, T}} P_0 \exp\left(\frac{L_\mathrm{v} \, {M} \, \left(T - T_\mathrm{v}\right)}{{R} \, T \, T_\mathrm{v}}\right), \end{equation} \tag{ 23 }$$

where P₀ is the ambient pressure, $L_\mathrm{v}$ the latent heat of vaporisation, M the molar mass, R the universal gas constant, T the temperature and $T_\mathrm{v} = 3315~\textrm{K}$ the evaporation temperature. The results obtained from the present numerical simulations using the heat source model introduced for Case 4 are compared with those reported by Cunningham et al [50] in figure 5, which show a reasonable agreement. The maximum absolute deviations between the present numerical predictions and experimental data for the melt-pool aspect ratio is less than 12%.

Zoom In Zoom Out Reset image size

Particular attention was also paid to the effect of spurious currents on numerical predictions of molten metal flow in melt pools with free surface deformations. Spurious currents are an unavoidable numerical artefact generating unphysical parasitic velocity fields that do not vanish with grid refinement when using the continuum surface force technique for modelling surface tension effects in the VOF method. This problem has been thoroughly investigated by Mukherjee et al [54] and the results of their study show that the Reynolds number based on the time-averaged maximum spurious velocity is $\mathcal{O}(10)$. The results of the present simulations show that the flow Reynolds number is $\mathcal{O}(10^{3})$ that is at least two orders of magnitude larger than that induced by spurious currents. Furthermore, to suppress fluid velocities in the solid regions a momentum sink term is defined based on the enthalpy porosity technique that further weakens the spurious currents in solid regions. Therefore, the influence of spurious currents on numerical predictions is considered to be negligible in the present study.

When accounting for surface deformations, the volume-of-fluid (VOF) method developed by Hirt and Nichols [14], based on a Eulerian formulation, is the most common method for modelling the melt pool behaviour. In this diffuse boundary method, the interfacial forces and the energy fluxes applied on the melt-pool surface are treated as volumetric source terms in the surface region, instead of imposing them as boundary conditions. In this approach, however, the fact that surface deformations lead to temporal and spatial variations of the free surface boundary conditions, as remarked by Meng et al [15] and Wu et al [16], is often neglected. The results reported by Choo et al [17] suggest that variations in power-density distribution and changes in free-surface profile can affect molten metal flow in melt pools and its stability. In this section, the influence of applying each of the five different approaches to model the laser heat source (described as Cases 1–5 in section 2.2.2) on the melt pool behaviour is discussed. In these simulations, no enhancement factor, as introduced in the previous section, was employed since it is ad hoc and has little justification in physical reality.

The energy flux distribution over the melt-pool surface determines the spatial temperature distribution over the melt-pool surface and its temporal variations. Due to the temporal changes in the melt-pool surface shape, the total energy absorbed by the material in Case 1 varies between $491.4 \, \textrm{W}$ and $538.1 \, \textrm{W}$, with a median value of $515.7 \, \textrm{W}$, which differs from the total energy supplied by the laser $\eta \, Q = 500.5 \, \textrm{W}$. This issue is resolved by utilising a dynamic adjustment coefficient in Case 2, while the relative energy flux distribution remains unchanged. However, the results reported by Courtois et  al [55] and Bergström et  al [56] showed that the energy flux distribution varies with surface deformation during a laser melting process. When surface deformations are too small to cause multiple reflections, redistributing the energy flux over the free surface results in a better input energy conservation as obtained in Case 3. The total energy absorbed by the material in Case 3 ranges between $482.6 \, \textrm{W}$ and $532.3 \, \textrm{W}$, with a median value of $503.8 \, \textrm{W}$. By further introducing an adjustment coefficient (equation (17)), the absorbed energy can be made to exactly match the supplied laser power ηQ. Finally, for the most simple approach without any surface deformations (i.e., Case 5), the total absorbed energy exactly matches the supplied laser power and remains unchanged in time.

The variations of temperature distribution over the melt-pool surface over time determine the thermocapillary driven flow pattern and melt-pool shape, as shown in figure 6. An unsteady, asymmetric, outwardly directed fluid flow emanating from the pool centre is found for all 5 cases. Flow accelerates towards the pool rim, where it meets inwardly directed fluid flow from the rim, due to the change of sign of $\left(\partial \sigma / \partial T\right)$ at a certain pool temperature (see figure 2). Close to the region where these two flows meet, velocity magnitudes are high due to the large thermocapillary stresses generated by the steep temperature gradients, in contrast to the low flow velocities in the central region of the pool surface. Interactions between these two opposing flows, and their interactions with the fusion boundary, cause the fluid flow inside the melt pool to be asymmetric and unstable [20], leading to a distorted melt-pool shape.

Zoom In Zoom Out Reset image size

Melt-pool surface temperatures are roughly 7–16% higher when surface deformations are taken into consideration (Cases 1–4), compared to flat surface simulations (Case 5). Temperature gradients over the melt-pool surface are also on average larger in Cases 1–4 (with a deformable surface) and particularly in Cases 2–4 (with heat source adjustment and/or redistribution), compared to Case 5 (with a non-deformable surface). Consequently, local fluid velocities caused by thermocapillary stresses are almost 25% higher for Cases 1–4, compared to those of Case 5, and are of the order of $0.6 \, \textrm{ms}^{-1}$, in reasonable agreement with experimental measurements performed by Aucott et  al [8] and estimated values from the scaling analyses also reported by Oreper and Szekely [2], Rivas and Ostrach [57] and Chakraborty and Chakraborty [58].

Free surface deformations have a destabilising effect on the melt pool due to the augmentation of temperature gradients over the melt-pool surface, leading to higher thermocapillary stresses [59]. Additionally, the stagnation region, where the sign of $\left(\partial \sigma / \partial T\right)$ and thus thermocapillary flow direction change, is sensitive to small spatial disturbances and further enhances melt pool instabilities. Furthermore, variations in magnitude and direction of velocities can cause rotational and pulsating fluid motions, leading to cross-cellular flow patterns with a stochastic behaviour in the melt pool [12, 60].

All such flow instabilities reinforce unsteady energy transport from the melt pool to the surrounding solid material, resulting in continuous melting and re-solidification of the material close to the solid-liquid boundary. This complex interplay leads to the melt-pool surface width for Case 5 to be 23% smaller than that for Case 1, and 10% smaller than for Cases 2–4. The higher amount of heat absorbed by the deformed free-surface in Case 1, and the enhanced convective heat transfer in Cases 1–4, are the main causes for the observed pool size differences.

To investigate the rotational fluid motion, the spatially-averaged angular momentum over the melt-pool surface ($\mathscr{L}$) about the z-axis (i.e. the axis of rotation) with respect to the origin $\mathscr{O}(0, 0, 0)$ is calculated as follows:

$$\begin{equation} \mathscr{L}\left(t\right) = \frac{1}{A}\iint \limits_{\mathrm{free\ surface}} \hat{z} \cdot \left(\vec{r} \times \vec{V}\right) \mathop{}\mathrm{d}A, \end{equation} \tag{ 24 }$$

where $\hat{z}$ is a unit vector in the z-direction, $\vec{r}$ the position vector, $\vec{V}$ velocity vector and A the area of gas-metal interface. Figure 7 shows the temporal variations of angular momentum over the melt-pool surface as a function of time. Positive values of the angular momentum in figure 7 show a clockwise flow rotation, and vice versa. Kidess et al [20] applied a similar approach to investigate rotational fluid motion over a flat non-deformable melt-pool surface. Continuous fluctuations in the sign of angular momentum indicate an oscillatory rotational fluid motion in the pool that is indeed self-excited. In Cases 1–4, flow pulsations start to take place after roughly $1 \, \textrm{s}$ and continue thereafter. This is not however valid for Case 5, where flow pulsations start to occur already after about $0.1 \, \textrm{s}$, and a clockwise rotation is established after $1.8 \, \textrm{s}$, which lasts for roughly two seconds. These instabilities in the flow field are attributed to a large extent to the interactions between the two opposing flows meeting at the melt-pool surface [13, 20], and to a lesser extent the effects of hydrothermal waves [61, 62].

Zoom In Zoom Out Reset image size

Free surface deformations evolve rapidly because of the thermocapillary stresses acting on the melt-pool surface and are shown in figure 8 at three different time instances. Free surface deformations are smaller and less intermittent in Case 1 compared with Cases 2–4. This is due to the smaller thermocapillary stresses and the wider stagnation region, which makes the surface flow field less sensitive to spatial disturbances. Additionally, the Capillary number ($\mathrm{Ca} = \mu \lVert\vec{V}\rVert / \sigma$), which represents the ratio between viscous and surface-tension forces and which is $\mathcal{O}(10^{-3})$ for all the cases studied in the present work, appears to be larger in Cases 2–4 ($\mathrm{Ca} \approx 2.5\times 10^{-3}$) compared with Case 1 ($\mathrm{Ca} \approx 1.5\times 10^{-3}$), particularly in the central region of the melt-pool surface. Redistributing the energy flux based on the free surface profile (Cases 3 and 4) disturbs the temperature field locally and thus the velocity distribution over the melt-pool surface.

Zoom In Zoom Out Reset image size

Fluid flow may influence the energy transport and associated phase change significantly during laser melting [7], which can be assessed through the Péclet number that is the ratio between advective and diffusive heat transport, defined as follows:

$$\begin{equation} \mathrm{Pe} = \frac{\mathscr{D} \lVert\vec{V}\rVert \rho c_\mathrm{p}}{k}, \end{equation} \tag{ 25 }$$

where $\mathscr{D}$ is a characteristic length scale, here chosen to be $r_\mathrm{b}$. The value of the Péclet number is much larger than one ($\mathrm{Pe} = \mathcal{O}(10^2)$), which indicates a large contribution of advection to the total energy transport, compared with diffusion.

The predicted melt-pool depth is shown in figure 9 at two different time instances. When surface deformations are taken into account (Cases 1–4), the predicted melt-pool shape is shallow, with its maximum depth located in the central region and its largest width at the surface. When the surface is assumed to be non-deformable (Case 5), the depth profile of the melt-pool resembles a doughnut-shaped, with its maximum width located below the gas-metal interface. The predicted fusion boundaries of all the cases studied in the present work are asymmetric and unstable. Flow instabilities in the melt pool increase when surface deformations are taken into account, which enhance convection in the pool resulting in a smooth melt-pool depth profile. Variations in the melt-pool shape for Case 5 are less conspicuous compared to those for Cases 1–4 because of the smaller fluctuations in the fluid flow field.

Zoom In Zoom Out Reset image size

To understand the three-dimensional oscillatory flow behaviour, three monitoring points were selected over the melt-pool surface in different non-azimuthal directions $p_1\left(x^*, y^*\right) = \left(5/7,0\right)$, $p_2\left(x^*, y^*\right) = \left(5/7,5/7\right)$ and $p_3\left(x^*, y^*\right) = \left(0,5/7\right)$ in addition to a point $p_4\left(x^*, y^*\right) = \left(0,0\right)$ on the melt-pool surface, in which x^* and y^* are non-dimensionalized with the laser beam radius $r_\mathrm{b}$. Temperature signals, recorded from the monitoring point p₄, and the corresponding 'fast Fourier transform' (FFT) [63] frequency spectra are shown in figure 10. In the presence of surface deformations (Cases 1–4), self-excited temperature fluctuations, initiated by numerical noise, grow and reach a quasi-steady state after about $1 \, \textrm{s}$. Without surface deformations (Case 5), the amplitudes of temperature fluctuations remain very small. In Case 1, where surface deformations are taken into account, the temperature fluctuations with large amplitudes have a fundamental frequency of $f_\mathrm{0} \approx 22~\textrm{Hz}$, with harmonics at $f_1 \approx 2 f_\mathrm{0}$ and $f_2 \approx 3 f_\mathrm{0}$. Irregular patterns appear in the spectrum of temperature fluctuations for Case 2–4, showing a highly unsteady behaviour that results from the enhancement of thermocapillary stresses and the variations in energy flux distribution.

Zoom In Zoom Out Reset image size

Temperature signals recorded from monitoring points p₁, p₂ and p₃ in the time interval of 2–$3 \, \textrm{s}$ are shown in figure 11. The data for Case 1 indicate that the thermal and fluid flow fields are dominated by a pulsating behaviour. This is also valid for Case 4 up to roughly $2.5 \, \textrm{s}$; however, after $2.5 \, \textrm{s}$ modulation of temperature fluctuations takes place, which is probably due to changes in temperature gradients, melt pool size and the complex interactions between vortices generated inside the melt pool. Similar behaviour was observed for Cases 2 and 3 (not shown here). The amplitudes of temperature fluctuations appear to be smaller for Case 5 compared to the other cases and show an irregular behaviour, revealing complex unsteady flow in the melt pool.

Zoom In Zoom Out Reset image size

In this section, the effects of employing temperature-dependent properties on numerical predictions of thermal and flow fields as well as the melt-pool shape are investigated using two different heat source models (Cases 1 and 4). Temperature-dependent properties for the metallic alloy considered in the present study (S705) are given in table 4, where the values are estimated by analogy with the values for iron-based alloys [64]. Temperature distribution over the melt-pool surface and free-surface flow after $5 \, \textrm{s}$ of heating are shown in figure 12 for both temperature-dependent and temperature-independent thermophysical properties. Melt-pool surface temperatures predicted using temperature-independent properties are about 1–6% lower than those predicted using temperature-dependent properties. This is mainly attributed to the changes in the thermal diffusivity $\left(\alpha = k / \left(\rho c_\mathrm{p}\right)\right)$ of the molten material. The thermal conductivity of the molten material changes from $21.8\,\textrm{W}\,\textrm{m}^{-1}\,\textrm{K}^{-1}$ at $T = 1620 \, \textrm{K}$ to $33.4\,\textrm{W}\,\textrm{m}^{-1}\,\textrm{K}^{-1}$ at $T = 2700 \, \textrm{K}$. Hence, the average thermal diffusivity of the molten material obtained from a temperature-dependent model is about 10% lower than that estimated using temperature-independent properties. The fluid velocities are roughly 15–32% lower when temperature-independent properties are employed, compared to those predicted using temperature-dependent properties, which is due to the decrease in the viscosity of the molten material at elevated temperatures. With an increase in temperature, the thermal conductivity of the molten metal increases and the viscosity of the molten metal decreases, resulting in a reduction of momentum diffusivity and enhancement of thermal diffusivity and thus reduction of the Prandtl number $\left(\mathrm{Pr} = c_\mathrm{p} \, \mu / k\right)$. Although the heat source models affect the thermal and flow fields in the melt pool, it is found that the numerical predictions are less sensitive to the heat source models when temperature-dependent properties are employed for the cases studied in the present work. However, it should be noted that for the cases where surface deformations are larger, the effects of heat source adjustment on numerical predictions become critical, as discussed in section 4.2.

Zoom In Zoom Out Reset image size

Table 4. Temperature-dependent thermophysical properties of the Fe–S alloy used in the present study. Values are estimated by analogy with the values for iron-based alloys [64].

Property	Fe-S alloy	Unit
Density ρ	8100	kg m−3
Specific heat capacity $c_\mathrm{p}$	627 (solid phase)	J kg−1 K−1
	723.14 (liquid phase)
Thermal conductivity k	8.8521+0.0114·T (solid phase)	W m−1 K−1
	4.5102+0.0107·T (liquid phase)
Viscosity µ	${0.0659} - {7\times 10^{-5}}\cdot T + {3\times 10^{-8}}\cdot T^2 - {3\times 10^{-12}}\cdot T^3$	kg m−1 s−1
Latent heat of fusion $L_\mathrm{f}$	250 800	J kg−1
Thermal expansion coefficient β	2 × 10−6	K−1
Liquidus temperature $T_\mathrm{l}$	1620	K
Solidus temperature $T_\mathrm{s}$	1610	K

Changes in the material properties with temperature affect thermal and flow fields in the melt pool, resulting in changes in the predicted melt-pool shape, as shown in figure 13. When temperature-dependent properties are employed, the melt-pool depth is significantly (about 36–74%) larger than that predicted using temperature-independent properties.

Zoom In Zoom Out Reset image size

Temperature signals recorded from the monitoring point p₄ and the corresponding frequency spectra are shown in figure 14. Temperature signals received from the monitoring point p₄ are almost in the same range and vary between $2320 \, \textrm{K}$ and $2820 \, \textrm{K}$. Employing temperature-dependent properties affects the frequency spectra of fluctuations, however the amplitude of fluctuations remains almost unaffected. The results presented in figure 14 indicate that adjusting the heat source dynamically during simulations can result in a decrease in the amplitude of temperature fluctuations, however its effect on the frequency spectrum is insignificant.

Zoom In Zoom Out Reset image size

In section 4.1 it was mentioned that, in many simulation studies in literature, the thermal conductivity and viscosity of the liquid metal were artificially increased by a so-called enhancement factor $\mathcal{F}$, in order to obtain better agreement between experimental and simulated post-solidification melt-pool shapes. Numerical studies carried out by De et al [65–67] showed that the values reported for the enhancement factor in the literature depend greatly on operating conditions and ranges from 2 to 100 (see for instance, [68–71]), however values between 2 and 10 are most often employed. The effects of the enhancement factor $\mathcal{F}$ on thermal and fluid flow fields in molten metal melt pools as well as the melt-pool shape can be found in [19, 43, 45], thus are not repeated here. Focusing on Case 4 in which free surface deformations are accounted for, the influence of $\mathcal{F}$ on the oscillatory flow behaviour is investigated. Temperatures recorded from the monitoring point p₄ and the temperature fluctuation spectra are shown in figure 15. Temperatures and the amplitudes of fluctuations reduce with increasing $\mathcal{F}$. The contribution of diffusion in total energy transport increases with $\mathcal{F}$, which results in a reduction of temperature gradients and therefore thermocapillary stresses generated over the melt-pool surface. The reduced thermocapillary stresses in addition to the enhanced viscosity of the molten metal lead to a reduction of fluid velocities that decrease convection in the melt pool further, which significantly affects the fluid flow structure. Increasing $\mathcal{F}$ to 2.5 results in relatively deeper melt pool with a corresponding reduction of the fundamental frequency of fluctuations. Using higher values of $\mathcal{F}$, the melt-pool shape approaches a spherical-cap shape and the frequency spectrum becomes rather uniform with small amplitude fluctuations.

Zoom In Zoom Out Reset image size