Comparative Numerical Analysis of Keyhole Shape and Penetration Depth in Laser Spot Welding of Aluminum with Power Wave Modulation

Abstract

Keyhole mode laser welding is a valuable technique for welding thick materials in industrial applications. However, its susceptibility to fluctuations and instabilities poses challenges, leading to defects that compromise weld quality. Observing the keyhole during laser welding is challenging due to bright process radiation, and existing observation methods are complex and expensive. This paper alternatively presents a novel numerical modeling approach for laser spot welding of aluminum through a modified mixture theory, a modified level-set (LS) method, and a thermal enthalpy porosity technique. The effects of laser parameters on keyhole penetration depth are investigated, with a focus on laser power, spot radius, frequency, and pulse wave modulation in pulsed wave (PW) versus continuous wave (CW) laser welding. PW laser welding involves the careful modulation of power waves, specifically adjusting the pulse width, pulse number, and pulse shapes. Results indicate a greater than 80 percent increase in the keyhole penetration depth with higher laser power, pulse width, and pulse number, as well as decreased spot radius. Keyhole instabilities are also more pronounced with higher pulse width/numbers and frequencies. Notably, the rectangular pulse shape demonstrates substantially deeper penetration compared to CW welding and other pulse shapes. This study enhances understanding of weld pool dynamics and provides insights into optimizing laser welding parameters to mitigate defects and improve weld quality.

1. Introduction

In recent years, high-power laser welding has been extensively utilized in diverse industrial sectors, such as automotive manufacturing, astronautics and aeronautics, pipelines, and reactor vessels [1,2], to name a few. High-power laser welding is a promising candidate for industrial use due to its potential to efficiently join thick plates of various metals, including aluminum and steel. It enables single-pass welding of thick plate metals while considerably increasing the welding efficiency [3]. There are two laser welding modes: conduction and keyhole. In conduction mode, low energy power density deposition through shaped beam distribution leads to material fusion, but with weak penetration. Conversely, the keyhole mode benefits from high energy power density (above 10⁶ W cm⁻²), raising the material’s temperature above its vaporization point. This leads to metal evaporation from the material surface, during which a keyhole is formed inside the material due to the vaporization pressure produced, also known as recoil pressure; the result is a deep weld penetration, ranging from millimeters and extending to centimeters [4,5]. Different laser types, including continuous wave (CW), pulsed wave (PW), and modulated wave (MW) lasers, can be employed for welding, depending on the material being used. Aluminum and its alloys possess desirable properties such as oxidation and corrosion resistance, thermal and electrical conductivity, high reflectivity, low density, and strength. PW lasers effectively overcome the reflectivity of aluminum by delivering high energy density during pulse initiation while controlling the average power to minimize the overall heat. PW lasers are thus suitable options for welding aluminum and its alloys [6,7,8]. Laser welding technology offers benefits such as high production rates, low energy input, faster welding speeds, smaller heat-affected zones, and narrower weld beads [9]. However, monitoring and controlling defects resulting from improper laser welding is crucial in order to prevent environmental and economic consequences [9,10]. Laser welding is associated with several defects, including humping [11], porosities [12], spatters [13], and other imperfections induced by a small beam diameter. Additionally, geometric defects such as a lack of fusion [14], sagged welds [15], and concave root surfaces [16] can arise from imprecise beam positioning or pore formation during deep penetration welding as a result of high cooling rates [9,10,17]. These defects considerably undermine the reliability of manufactured products [18]. According to the literature, the majority of these defects arise from the dynamic variations in the keyhole, keyhole collapse, and melt pool oscillations. In particular, porosities resulting from keyhole instabilities are prevalent in aluminum welding, highlighting the importance of controlling keyhole and melt pool fluctuations for welding process stability and quality [2,19]. Scientists use experimental observations and numerical simulations to monitor laser welding and keyhole dynamics. However, high-quality laser welds performed through experiments often rely on trial and error, which is time-consuming, costly, and not environmentally friendly. The challenge is amplified when welding materials with high thermal conductivity or volatile chemical elements such as aluminum and its alloys. Experimentation also fails to accurately identify an efficient window of welding parameters due to the severe fluctuations of the keyhole and molten pool during the welding process [20,21]. On the other hand, the time scale for keyhole formation and laser welding is in the order of one microsecond (10⁻⁶ s), making direct experimental observations difficult and invisible to the naked eye [22]. Researchers have also used high-speed cameras to observe keyhole dynamics, but are limited to capturing surface phenomena of the fusion zone and vapor plume above the keyhole, missing important vapor plume dynamics inside. However, understanding keyhole fluctuations is closely tied to the behavior of vapor plumes, as the dynamics of vapor plumes can influence the stability and morphology of the keyhole [23]. Alternatively, numerical simulations offer significant opportunities to model and comprehend the underlying mechanisms representing keyhole dynamics, which enables the prediction of defect formation [24,25]. To get a better grasp of the mechanisms responsible for defect formation, numerical models should take into account fluid flow in the fusion zone, and determine the dynamic keyhole shape, which is tightly coupled with phase change, heat transfer, and fluid flow [26]. One of the early studies in this context was conducted by Lee et al. [17], who proposed a model of heat and fluid flow to describe spot laser welding. Their study integrated the Volume of Fluid (VOF) method into their model, addressing vapor phase and molten pool surface displacement, while also investigating dynamic instabilities and keyhole oscillations. In another study, Medal et al. [27] developed a model for stationary laser welding and used a moving mesh method (Arbitrary Lagrangian–Eulerian) to track the free surface. Even though they considered recoil pressure and vaporization impacts at the free surface through specific boundary conditions, the presence of mesh distortion hindered the prediction of large displacements beyond the early stages of keyhole formation. Pang et al. [28] used the LS method to expand a model for predicting keyhole formation and studied keyhole instability mechanisms under different heat input conditions. Courtois et al. [10,26] presented a comprehensive heat and fluid flow model for laser welding, incorporating all three phases (liquid, solid, and vaporized metal) while employing the LS method to track the liquid/vapor interface. Their model effectively predicted keyhole oscillations and porosity formation, with acceptable agreements with experimental results. However, the significantly higher laser wavelength used in their research, compared to what is typically used in industrial welding, raises concerns about its applicability. Furthermore, the omission of the Marangoni effect, driven by temperature-induced surface tension gradients, might impact keyhole dynamics and the overall welding process, thereby forcing a further assessment of these factors. Matsunawa and Katayama [29,30] carried out a series of experimental investigations studying the impact of keyhole instability in deep penetration laser welding. Their findings corroborated the correlation between porosity formation and keyhole instabilities, which has not been examined in other welding fusions. The study disclosed an unstable keyhole in CW, exhibiting dynamic variations in size, shape, and depth. This led to bubble formation at the bottom part of the molten pool. An obstruction of the upward bubble flow by solidified metal resulted in keyhole-induced porosity. To control and address keyhole instability and its significant impact on defect formation, various methods have been proposed, such as pulse modulation [29], the use of side assisting gas [31], dual-beam welding [21], and beam oscillation [32]. Ke et al. [32] numerically studied keyhole-induced porosity in laser beam oscillating welding of 5A06 aluminum alloy. They evaluated the impact of different laser paths, including oscillating and non-oscillating paths, using a hybrid heat source model comprising an adaptive Gaussian profile with a double ellipsoid heat source. The results indicated that beam oscillation contributed to the formation of broader and more stable keyholes, a shallower and larger molten pool, and a more complex fluid flow, which minimized the likelihood of keyhole collapse, bubble formation, and porosity formation. Pulse modulation is another technique that can reduce the oscillation and instability of molten pools and keyholes. Matsunawa et al. [29] extensively investigated keyhole and molten pool behavior, revealing that keyhole stability is adversely affected by intense metal evaporation on the keyhole wall at its front position. This leads to porosity formation inside the weld the weld pool and in the final product. Their findings indicated that keyhole-induced porosity could be significantly suppressed by employing pulse modulation. By implementing proper duty cycles and frequencies in pulse modulation, the number of porosities can be reduced through the effective removal of holes benefiting from subsequent pulses and appropriate overlapping ratios.

It can be seen in the literature that the energy density of the laser plays a crucial role in keyhole formation and behaviors. To the best of our knowledge, the impact of various laser energy pulse shapes on keyhole mode laser welding has not yet been covered in the literature; also, current pulse shapes are mostly confined to rectangular pulse shapes. The present study aims to investigate the impact of various laser energy pulse shapes (trapezium, triangular, rectangular, and variant–rectangular) on the keyhole and molten pool shape in a 2D axisymmetric configuration case. The results of keyhole penetration depth and morphology for these pulse shapes are compared with those from continuous laser welding to analyze the impact of pulse modulation on keyhole formation. Moreover, aluminum was chosen as the base metal for the investigation, as laser welding on it presents challenges in terms of numerical convergencies, and only a very limited number of numerical works have explored laser welding on the metal. Furthermore, given that aluminum is much lighter than other metals such as steel, and given its high thermal conductivity, which makes it highly reflective, studying its weldability and improving its welding efficiency are a must, and greatly relevant to automotive manufacturing. Correspondingly, a novel numerical model equipped with modified techniques, such as modified mixture theory and the modified LS method, is utilized to investigate the phase transformations and coupled physics in this problem. The former defines the material properties in the entire domain, including the gas, solid, and liquid states, while the latter tracks the vapor/liquid interface inside the material due to the evaporation phenomena during laser welding. Additionally, the model developed in this paper is validated by comparing it with experimental results obtained by Qin et al. [33], shown in Section 4.1 of the Results and Discussion Section.

2. System Description and Material

The present paper presents a comparative numerical study of keyhole dynamics and the molten pool shape during laser spot welding of aluminum. A 2D axisymmetric model was utilized to represent the base metal. Its dimensions and configuration are depicted in Figure 1a. The choice of a 2D axisymmetric model for the simulation is rationalized based on the stationary spot laser welding assumption in this work. The laser beam position remains fixed within the coordinate system during the spot laser application, making the laser welding rotationally symmetric with respect to the vertical z-axis. This configuration thus represents a transversal cross-section of the sample for keyhole monitoring, as shown in Figure 1b. The plane used to monitor the keyhole morphology and depths in the Result and Discussion Section is also depicted in Figure 1c. It is worth mentioning that all the simulations, developments, models, and the investigated results of keyhole morphology and depth in the present paper were built in the graphical user interface of COMSOL Multiphysics 5.6. Using this software, the transition from a 2D axisymmetric model to a 3D model is accessible immediately. While COMSOL simplifies the process and offers 3D simulation results for our model, it is important to note that the software automatically imposes symmetry conditions for axisymmetric assumptions during calculations. This means it disregards variations and fluctuations occurring outside the plane of symmetry. Consequently, this excludes three-dimensional flow patterns or convective currents, which play a role in influencing heat transfer and material transport dynamics within the weld pool region. Therefore, the model could not be considered a fully comprehensive 3D representation and does not fully mimic industrial scenarios. Nonetheless, it serves as a crucial initial step towards 3D modeling, aiding in the validation of the numerical model.

The material used as the base metal was aluminum and its properties are presented in

Table 1. Thermophysical properties of aluminum [34,35,36].

Property	Symbol	Magnitude
Solidus temperature	$T_s$	847 (K)
Liquidus temperature	$T_l$	905 (K)
Vaporization temperature	$T_V$	2743 (K)
Solid density	${\rho }_s$	2700 (kg/m3)
Liquid density	${\rho }_l$	2385 (kg/m3)
Solid thermal conductivity	$k_s$	238 (W/m/K)
Liquid thermal conductivity	$k_l$	100 (W/m/K)
Liquid specific heat capacity	$C_{p,s}$	917 (J/kg/K)
Solid specific heat capacity	$C_{p,l}$	1080 (J/kg/K)
Latent heat of fusion	$L_m$	3.896 × 105 (J/kg)
Latent heat of vaporization	$L_V$	9.462 × 106 (J/kg)
Radiation emissivity	$\xi$	0.2
Convective heat transfer coefficient	h	20 (W/m2/K)
Thermal expansion coefficient	$\beta$	2.36 × 10−5 (1/K)
Dynamic viscosity	$\mu$	1.6 × 10−3 (Pa.s)
Surface tension coefficient	$\sigma$	0.95 × (1 + 0.13 × (1 − T/Tm))1.67 (N/m)
Surface tension coefficient with temperature	$\partial \sigma /\partial T$	−0.3 × 10−3 (N/m/K)

. Moreover, a series of tests were designed to investigate the impact of various laser and process parameters. All the investigated cases are presented in

Table 2. List of the investigated sample cases with related operating properties for variations in laser characteristics (LC) and for modulated wave (MW) welding.

Case No.	Laser Power	Pulse Width	Number of Pulses	Frequency of Laser	Period of Pulse	Pulse Shape	Spot Radius	Total on Time
LC1	6 kW	2 ms	1	100 Hz	0.01 s	Rectangular	300 µm	0.002 s
LC2	6 kW	2 ms	1	100 Hz	0.01 s	Rectangular	425 µm	0.002 s
LC3	6 kW	2 ms	1	100 Hz	0.01 s	Rectangular	525 µm	0.002 s
LC4	6 kW	2 ms	1	100 Hz	0.01 s	Rectangular	725 µm	0.002 s
Impact of spot radius
LC5	6 kW	1 ms	3	50 Hz	0.0066 s	Rectangular	300 µm	0.003 s
LC6	6 kW	1 ms	3	100 Hz	0.0033 s	Rectangular	300 µm	0.003 s
LC7	6 kW	1 ms	3	150 Hz	0.0022 s	Rectangular	300 µm	0.003 s
Impact of frequency
LC8	2 kW	3 ms	1	100 Hz	0.01 s	Rectangular	300 µm	0.003 s
LC9	4 kW	3 ms	1	100 Hz	0.01 s	Rectangular	300 µm	0.003 s
LC10	6 kW	3 ms	1	100 Hz	0.01 s	Rectangular	300 µm	0.003 s
Impact of laser power
MW1	6 kW	0.5 ms	1	100 Hz	0.01 s	Rectangular	300 µm	0.0005 s
MW2	6 kW	1 ms	1	100 Hz	0.01 s	Rectangular	300 µm	0.001 s
MW3	6 kW	2 ms	1	100 Hz	0.01 s	Rectangular	300 µm	0.002 s
MW4	6 kW	3 ms	1	100 Hz	0.01 s	Rectangular	300 µm	0.003 s
Impact of Pulse width
MW5	4 kW	0.5 ms	2	100 Hz	0.005 s	Rectangular	300 µm	0.001 s
MW6	4 kW	0.5 ms	6	100 Hz	0.0016 s	Rectangular	300 µm	0.003 s
MW7	4 kW	0.5 ms	10	100 Hz	0.001 s	Rectangular	300 µm	0.005 s
MW8	4 kW	0.5 ms	14	100 Hz	0.00071 s	Rectangular	300 µm	0.007 s
MW9	4 kW	0.5 ms	18	100 Hz	0.00055 s	Rectangular	300 µm	0.009 s
Impact of pulse number
CW	2 kW	10 ms	1	100 Hz	0.01 s	Continuous	300 µm	0.01 s
Impact of continuous laser welding
MW10	4 kW	5 ms	1	100 Hz	0.01 s	Rectangular	300 µm	0.005 s
MW11	4 kW	8 ms	1	100 Hz	0.01 s	Trapezium	300 µm	0.008 s
MW12	4 kW	10 ms	1	100 Hz	0.01 s	Triangle	300 µm	0.01 s
MW13	4 kW	8 ms	1	100 Hz	0.01 s	Trap.: t2	300 µm	0.008 s
MW14	1–3 kW	Variant	1	100 Hz	0.01 s	Var.-Rect.	300 µm	0.01 s
MW15	2–4 kW	Variant	1	100 Hz	0.01 s	Rect.-Tri.	300 µm	0.008 s
MW16	2–3 kW	Variant	1	100 Hz	0.01 s	Rect.-Trap.	300 µm	0.008 s
MW17	2–3 kW	Variant	1	100 Hz	0.01 s	Rect.-Rect.	300 µm	0.008 s

Impact of pulse shape (identical total 20 J laser power was used for cases MW10-17 and CW), Var.: variant, Rect.: Rectangular, Trap.: Trapezium, Tri. Triangle, t2: type 2.

. Both CW and PW laser welding were examined for a detailed comparison. The PW laser welding was initially studied under different laser characteristics (LC), with the impact of various spot radii, laser frequency, and laser power on the keyhole penetration depth and shape investigated. Then, an investigation was performed into MW laser welding to examine the effects of pulse modulation, specifically, through the power modulation of waves, on the keyhole penetration depth. Various parameters, such as the pulse number, pulse width, and pulse shapes, were looked at. The pulse shapes selected for this study included rectangular, trapezium, and triangular shapes, as well as multi-shape pulses such as rectangular–triangular, rectangular–trapezium, and variant–rectangular pulses. Results were compared with those of CW laser welding to determine the cases exhibiting the greatest stability and keyhole penetration depth. The equivalent average energy of the laser was kept constant for the pulse shape impact, ensuring that investigated cases were physically comparable under identical operational conditions and heat source characteristics. The different pulse shapes utilized in the present paper are shown in Figure 2.

3. Numerical Approach and Methods

3.1. Heat and Fluid Flow Model

Laser heating occurs when the energy absorbed by the base metal from the laser beam converts into thermal energy. Subsequently, the temperature of the material surface increases due to the absorbed energy of the laser by the material. When the temperature surpasses the melting and vaporization thresholds, the material undergoes consecutive phase transitions, including fusion and evaporation, leading to keyhole formation [37]. The dynamics of the fluids are modeled by adding the impacts of recoil pressure, surface tension, and hydrostatic and hydrodynamic pressures [9,37]. Recoil pressure acts to open the keyhole, while the other forces act to close it, with the surface tension effect dominating the shape of the vapor/liquid interface, particularly at the top of the melt pool. The surface tension effect is attributable to two factors, namely, the Marangoni effect and the curvature effect, caused by the tangential and normal forces exerted on the free surface of the fluid [9,38]. Among the methods developed to model the free surface during a finite element computation, the modified LS method, as a Eulerian approach, is used in this study to model the free surface of the fluid and the vapor/liquid interface, as well as calculate the tangential and normal components of the free surface. Additionally, a new method called modified mixture theory is used to deal with phase transitions and discontinuities at the interfaces in the finite element computation [34,39]. The present paper aims to use the above approaches to represent and simulate the keyhole formation and the molten pool while taking into account the surface tension effects, free surface motion, buoyancy forces, recoil pressure, evaporation impacts, and mass loss due to evaporation. The following assumptions were made in the numerical investigation:

• The flow of molten material inside the fusion zone was assumed Newtonian, incompressible, and laminar.
• The temperature-dependent effects on the thermophysical properties and absorption coefficients were neglected for the sake of simulation simplicity.
• A porous medium, saturated with the liquid molten metal, was assumed for the mushy zone [34].
• A Gaussian laser beam distribution was assumed for the heat source.
• The impact of natural convection was added using the Boussinesq approximation [40].
• Plasma and the Knudsen layer were not taken into account.
• Multiple reflections of the laser beam were neglected in this model.
• The vaporized material known as metallic vapor was considered an ideal gas and transparent to the incoming laser beam.
• The thermal enthalpy porosity technique was used to track the solid/liquid interface and adds the impacts of temperature-dependent phase transitions (melting and vaporization) on the specific heat capacity in the heat transfer model [34].

3.2. Governing Equations

This section presents the Navier–Stokes and conservation equations in a standard format applicable to 2D axisymmetric configurations. It is noted that under the 2D axisymmetric assumption, the gradient of parameters is independent of the rotational angle. The transport phenomena in all three phases are calculated by simultaneously solving the equations of mass conservation, momentum, and energy conservation (Equations (1)–(3)):

$$\nabla .\overrightarrow{u}=0$$

(1)

$$\rho \left(\frac{\partial \overrightarrow{u}}{\partial t}+\overrightarrow{u}.\left(\overrightarrow{\nabla }.\overrightarrow{u}\right)\right)=\overrightarrow{\nabla }.[-pI+\mu (\overrightarrow{\nabla }\overrightarrow{u}+{(\overrightarrow{\nabla }\overrightarrow{u})}^T)])+\rho \overrightarrow{g}+{\overrightarrow{F}}_{Vol}$$

(2)

$$\rho C_p\frac{\partial T}{\partial t}+\rho C_p\overrightarrow{u}.\nabla T=\nabla .\left(k\nabla T\right)$$

(3)

where T represents the temperature, $\rho$ denotes material density, $C_p$ is the specific heat capacity, k stands for the thermal conductivity, and $\overrightarrow{u}$ the velocity vector. Moreover, I denotes the identity matrix, p represents the pressure, $\mu$ is the dynamic viscosity, and ${\left(\overrightarrow{\nabla }\overrightarrow{u}\right)}^T$ indicates the transpose of the gradient of the velocity vector $\overrightarrow{u}$. $\rho \overrightarrow{g}$ represents the effect of gravity, ${\overrightarrow{F}}_{Vol}$ is the total body force which comprises impacts of surface tension (${\overrightarrow{F}}_{st}$) and buoyancy effect (${\overrightarrow{F}}_{Buoyancy}$) through a Boussinesq approximation and Darcy damping force (${\overrightarrow{F}}_{Darcy}$), which are further described later in this work. The present paper uses various methods to model the melting/solidification and evaporation phase changes. The melting/solidification is modeled using the thermal enthalpy porosity technique [34,41], while the evaporation is modeled using the conservative modified LS method [42]. Moreover, modified mixture theory is also utilized to apply the mixture effects easing the finite element computation at the interface, especially for the elements combined with more than one phase [34,39]. Together with the LS method, modified mixture theory helps to deal with the discontinuities at the interfaces during the transition from one phase to another. It is noted that all these methods and the mixture properties were added to the heat transfer, laminar flow, and level-set interfaces of COMSOL individually. Then, to apply the coupling effects of these methods and interfaces to work simultaneously in favor of the calculation, a specific coupling interface of COMSOL Multiphysics named the Multiphysics interface was used. Correspondingly, the laminar flow, mixture properties, and level-set were coupled using the two-phase flow interface while the laminar flow and the heat transfer interfaces were coupled using the nonisothermal flow of the Multiphysics interface. The two coupling effects were then coupled under the Multiphysics interface of COMSOL. The methods used in this paper are given in detail in the subsequent sections.

3.2.1. Modified Mixture Theory

Some previous studies applied modified mixture theory for the solid and liquid phases [39], due to the coexistence of the solid and liquid phases in a single cell for some regions within the mushy zone. The method uses a mixture of material properties derived from both solid and liquid states for density, specific heat, conductivity, and dynamic viscosity, thereby referred to as a mixture theory. However, the present study proposes that the so-called mixture theory can be further developed and extended for application to all three phases in the same fashion. Domain 2 in Figure 1 contains all three phases following the occurrence of melting and vaporization when the keyhole is formed. Therefore, to account for the impact of each phase in the case of the coexistence of all three phases in a single cell, a mixture of properties derived from all three states (solid, liquid, gas) is also used for the material properties to help solve the governing equations.

Conduction Mode

In the initial stages of the welding process, the primary welding mode is the conduction mode, which consists of the solid and liquid phases. Here, the temperature is not high enough to induce evaporation, but rather only exceeds the melting temperature to initiate melting. Therefore, only the liquid and solid phases are available in this mode. To represent the material properties in this mode, a mixture of temperature-dependent properties of solid and liquid phases of the material is defined, and the domain is considered as a bulk mixture of solid/liquid phases [5]. Equations (4)–(7) are used to describe this mode:

$${\rho }_{sl}=V_s{\rho }_s+{{\rho }_{l}V}_l$$

(4)

$$k_{sl}=V_s{k}_s+{k_{l}V}_l$$

(5)

$${\mu }_{sl}=V_s{\mu }_s+{{\mu }_{l}V}_l$$

(6)

$${Cp}_{sl}=\frac{{\rho }_s}{{\rho }_{sl}}{V}_s{Cp}_s+\frac{{\rho }_l}{{\rho }_{sl}}{V}_l{Cp}_l$$

(7)

where $\rho$, $k$, $\mu$, $Cp$ are the density, thermal conductivity, dynamic viscosity, and specific heat capacity of the aluminum, respectively. Moreover, V is specified as the volume fraction of the material, which is distinguished by using s and l subscripts for the solid and liquid phases.

Transition and Keyhole Mode

The more the heat flux of the laser is irradiated on the surface of the aluminum, the more the temperature is increased until the temperature reaches the boiling temperature of the material. The aluminum therefore transitions from conduction to keyhole mode by exceeding the boiling temperature. At this stage, Domain 2 comprises solid, liquid, and gas phases of aluminum because a portion of the aluminum is vaporized and forms the keyhole [5]. Hence, to account for the material properties in this mode and solve the Navier–Stokes equations, a mixture of solid, liquid, and gas temperature-dependent properties of the material is defined. In this approach, each cell in the domain is treated as a liquid with properties of a mixture of solid, liquid, and gas phases.

This section presents the equations used to define the mixture properties of the computational domain. The notations g, l, and s in the following equations are used to specify the gas, liquid, and solid phases, where ${\rho }_{sl}$, and ${\rho }_{slg}$ denote the densities of the bulk mixture of solid/liquid and solid/liquid/gas, respectively. Equations (8)–(10) define the mixture properties for this mode:

$${\rho }_{slg}={\rho }_{sl}{V}_{f,2}+{\rho }_g{V}_{f,1}$$

(8)

$$k_{slg}=k_{sl}{V}_{f,2}+k_g{V}_{f,1}$$

(9)

$${Cp}_{slg}=\frac{{\rho }_{sl}}{{\rho }_{slg}}{V}_{f,2}{Cp}_{Al,eff}+\frac{{\rho }_g}{{\rho }_{slg}}{V}_{f,1}{Cp}_g$$

(10)

where $V_{f,1}$ and $V_{f,2}$ are respectively the volume fraction of gas (Domain 1) and the volume fraction of solid/liquid bulk (Domain 2). ${\rho }_{slg}$, $k_{slg}$, and ${Cp}_{slg}$ are the density, thermal conductivity, and specific heat capacity of the bulk mixture of gas/liquid/gas in Domain 2, respectively.

3.2.2. Tracking the Solid/Liquid Interface

Thermal Enthalpy Porosity Technique

To comprehend the definition of ${Cp}_{Al,eff}$, an enthalpy porosity technique commonly used to account for the phase change transition in fixed-grid techniques has been defined [34,41]. The technique defines an equivalent specific heat capacity ${Cp}_{Al,eff}$ to add the latent heat effects of fusion and vaporization in which the temperature interval for melting and vaporization are based on Equations (12) and (13), respectively. Taking two Gauss functions of $D_m$ and $D_V$ around the melting temperature ($T_m$) and vaporization temperature ($T_V$) and multiplying them by latent heats of fusion ($L_m)$ and vaporization ($L_V)$, the equivalent specific heat capacity is defined by adding the specific heat of solid/liquid mixture ${Cp}_{sl}$ to the equation, as presented in Equations (11)–(13). The smoothing interval for the fusion (${dT}_m$) is set to 29 K, based on the solidus and liquidus temperatures of the material ($\frac{T_l-T_s}{2}$) [43]. In this study, the smoothing interval of vaporization (${dT}_V$) is set to 50 K, as proposed by Tomashchuk et al. [44]:

$${Cp}_{eff}={Cp}_{sl}+L_m{D}_m+L_V{D}_V$$

(11)

$$(\mathrm{a})D_m=\frac{{exp}^{\left(\frac{-{\left({T-T}_m\right)}^2}{{{dT}_m}^2}\right)}}{\sqrt{\pi {{dT}_m}^2}}(\mathrm{b}){T_m-dT}_m\to {T_m+dT}_m$$

(12)

$$(\mathrm{a})D_V=\frac{{exp}^{\left(\frac{-{\left({T-T}_V\right)}^2}{{{dT}_V}^2}\right)}}{\sqrt{\pi {{dT}_V}^2}}\left(\mathrm{b}\right){T_V-dT}_V\to {T_V+dT}_V$$

(13)

The enthalpy porosity technique also helps to identify and localize the solid/liquid interface during fusion caused by laser welding. In this technique, the melt interface is treated implicitly using a quantity known as the volume fraction of liquid ($V_l$), which is presented in Equation (14). In fact, the partially solidified region known as the mushy zone is treated as a porous medium. For each of the cells, the liquid fraction represents the porosity inside, in which the porosity and the liquid velocity are zero for the fully-solidified regions. The flow in the mushy zone follows the Darcy law, which induces a Darcy damping force. This force is governed by a frictional dissipation within the mushy zone based on the Carman–Kozeny equation, Equation (17), for flow through a porous medium. Consequently, the more the porosity in the mushy zone diminishes, the more the velocity and permeability decrease, and the liquid velocity becomes zero for the fully-solidified regions [34,45]. The Darcy damping force can be defined in the momentum equation, Equation (2), using a source term specified in Equations (16) and (17), where b is used to avoid division by zero, d is a constant proportional to the dendrite size, which is considered constant and set to 10⁻² cm [46], and $V_l$ and $V_s$ are volume fractions of liquid and solid.

$$V_l=\left\{\begin{array}{c}1,T>T_l\\ \frac{T-T_s}{T_l-T_s},T_s\le T\le T_l\\ 0,T<T_s\end{array}\right.$$

(14)

$$V_S=1-V_l$$

(15)

$$F_{DarcyDamping}=-{\mu }_{l}K\overrightarrow{V}$$

(16)

$$K=\frac{\frac{180}{d^2}{\left(1-V_l\right)}^2)}{{V_l}^3+b}$$

(17)

3.2.3. Tracking the Vapor/Liquid Interface

Modified Level-Set Method

Due to the numerous driving forces and temperature-dependent material properties acting on it, the liquid/vapor interface is the most challenging phenomenon to tackle in keyhole mode laser welding. To address this challenge, this study proposes a modified format of the conservative LS method, which integrates the Volume of Fluid (VOF) and the narrow band LS method while incorporating a gas dynamic source term [42]. This source term considers the effects of evaporation and the boiling effect induced by mass loss and recoil pressure on the interface. The standard and general format of the transport equation of the LS method is presented in Equation (18):

$$\frac{\partial \varphi }{\partial t}+\overrightarrow{u}.\nabla \varphi +{\gamma }_{ls}\nabla .(\varphi \left(1-\varphi \right)\frac{\nabla \varphi }{\left|\nabla \varphi \right|}-{ϵ}_{ls}\nabla \varphi )=0$$

(18)

where ${\gamma }_{ls}$ is the reinitialization parameter related to the flow velocity and ${ϵ}_{ls}$ the parameter that controls the interface thickness. Furthermore, $\varphi$ is the level-set variable defined using a Heaviside function [42]. The LS variable varies smoothly between 0 and 1 within the interface layer and is set to 0.5 for the vapor/liquid interface, as shown in Equation (19). The definition of this variable throughout all elements of the computational domain, along with its transport using fluid flow calculations, helps to track the vapor/liquid interface and distinguishing between gaseous and condensed phases.

$$\varphi \left(x,y,t\right)=\left\{\begin{array}{c}0,{T>T}_V,y<{-ϵ}_{ls}\\ 0.5,T_l<T<T_V,y=0\\ 1,T<T_l,y>{ϵ}_{ls}\end{array}\right.$$

(19)

Using the level-set variable $\varphi$, the normal and tangential vectors on the interface are defined in Equations (20) and (21) [34,39].

$$\overrightarrow{n}=\frac{\nabla \varphi }{\left|\nabla \varphi \right|}$$

(20)

$$\overrightarrow{k}=\nabla .\frac{\nabla \varphi }{\left|\nabla \varphi \right|}$$

(21)

Moreover, a delta function of the $\varphi$ variable is defined in Equation (22) to facilitate smooth transitions and phase traversal at the interfaces during the finite element computation [34]. This parameter definition aids in minimizing discontinuities, particularly for temperature-dependent material properties and forces that should only apply at the interface, such as the surface tension, recoil pressure, and the laser heat source. By multiplying these forces along with any other terms with the delta function, they are confined to the vapor/liquid interface during calculations, ensuring they possess a non-zero value solely at that interface.

$$\delta =6\left|\varphi (1-\varphi )\right|\left|\nabla \varphi \right|$$

(22)

3.2.4. Definition of Source Terms and Driving Forces on the Interface

Recoil Pressure, Mass Loss, and Evaporative Source Term of Heat Flux

As mentioned before, when the laser beam irradiates on the material surface, a keyhole is formed as the temperature exceeds the vaporization point, leading to the formation of a vapor pressure known as recoil pressure. The mass loss rate due to evaporation is derived by writing Navier–Stokes, continuity, and energy equations across the Knudsen layer. The mass loss rate is then expressed by the Hertz–Langmuir relation [43], as presented in Equation (23):

$${\dot{m}}_{H-L}=\sqrt{\frac{M}{2\pi R}}\frac{P_{sat}\left(T\right)}{\sqrt{T}}(1-{\beta }_r)$$

(23)

where ${\beta }_r$ represents the retro-diffusion coefficient, while R and M denote the universal gas constant and the molar mass of the evaporated particles, respectively. Moreover, $P_{sat}$ is the saturated vapor pressure expressed based on the Clausius–Clapeyron law as presented in Equation (24) [43,47], where $P_{atm}$ is defined as the atmospheric pressure and $T_v$ the evaporation temperature at the atmospheric pressure.

$$P_{Sat}=P_{atm}\mathit{exp}\left[\frac{M{L}_v}{R{T}_v}\left(1-\frac{T_v}{T}\right)\right]$$

(24)

The literature commonly employs empirical formulas to describe the recoil pressure. Mayi et al. [43] introduced a formula for recoil pressure as $P_{rec}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.(1+{\beta }_r)P_{Sat}\left(T\right)$ while Lee et al. [17] and Geiger et al. [48] proposed that the recoil pressure can be determined as $P_{rec}\cong 0.54P_{atm}\exp (∆H_{LV}T-T_{LV}/RT{T}_{LV})$, where $∆H_{LV}$ represents the evaporation enthalpy. However, a more tangible approach is to add the impact of the recoil pressure into the continuity equation by introducing a source term, as presented by Zhang et al. [34] and Courtois et al. [26]. Due to the differences between the gaseous and liquid densities (${\rho }_l-{\rho }_v$) in the vicinity of the vapor/liquid interface, the incompressibility of the fluid phases cannot be satisfied on two opposite sides of the interface. Hence, a source term is added to the continuity equation using the delta function of the LS variable defined in Section 3.2.3 to smooth the transition between phases for density variations and also add the mass loss impact and the recoil pressure. Correspondingly, the source term will be non-zero only on the interface, as presented in Equation (25) for the modified continuity equation, and the incompressibility of phases is satisfied for all the regions far from the interface [26,34]:

$$\nabla .\overrightarrow{u}=\delta \left(\varphi \right){\dot{m}}_{H-L}\left(\frac{{\rho }_l-{\rho }_v}{{\rho }^2}\right)$$

(25)

Moreover, another source term is also added to the transport equations of the LS variable to enhance the impact of evaporation phenomena induced by mass loss only on the interface. The source term will smooth the transport of the LS variable in the finite element computation on both sides of the interface using the delta function of the LS variable, volume fractions of gaseous and liquid phases, and their corresponding densities. The modified form of the transport equations of the LS variable is presented in Equation (26).

$$\frac{\partial \varphi }{\partial t}+\overrightarrow{u}.\nabla \varphi -\delta \left(\varphi \right){\dot{m}}_{H-L}\left(\frac{V_{f,1}}{{\rho }_v}+\frac{V_{f,2}}{{\rho }_l}\right)+{\gamma }_{ls}\nabla .(\varphi \left(1-\varphi \right)\frac{\nabla \varphi }{\left|\nabla \varphi \right|}-{ϵ}_{ls}\nabla \varphi )=0$$

(26)

where $V_{f,1}$ and $V_{f,2}$ are the volume fraction of gas (Domain 1) and the volume fraction of solid/liquid bulk (Domain 2), respectively.

3.2.5. Definition of the Surface Tension Impact and Boussinesq Approximation

Given the role of the delta function of the LS variable in smoothing transitions between phases and the application of the driving forces only on the interface, the surface tension impact and the buoyancy force are added to the momentum equation as body forces. Equations (27) and (28) are presented for the surface tension force (curvature and Marangoni effect) and the buoyancy force, respectively. Equation (28) considers the buoyancy effect but only within the molten aluminum using the term $\varphi$, and the effect is neglected in the gas phase due to intense vapor velocity. Then, the final format of the Navier–Stokes equations is presented in Equation (29). It should be noted that the buoyancy force is only applied within the molten pool utilizing the LS variable $\varphi$, and its impact is negligible in the gas phase [10,26].

$${\overrightarrow{F}}_{st}=(\gamma .nk-{\nabla }_s\gamma .t)\delta \left(\varphi \right)$$

(27)

$${\overrightarrow{F}}_{Buoyancy=}-{\rho }_l{\beta }_l(T-T_{melting})\overrightarrow{g}\varphi$$

(28)

$$\rho \left(\frac{\partial \overrightarrow{u}}{\partial t}+\overrightarrow{u}.\left(\overrightarrow{\nabla }.\overrightarrow{u}\right)\right)=\overrightarrow{\nabla }.[-pI+\mu (\overrightarrow{\nabla }\overrightarrow{u}+{\left(\overrightarrow{\nabla }\overrightarrow{u}\right)}^T\left)\right])+\rho \overrightarrow{g}-{\rho }_l{\beta }_l(T-T_{melting})\overrightarrow{g}\varphi +-{\mu }_{l}K\overrightarrow{V}+(\gamma .nk-{\nabla }_s\gamma .t)\delta \left(\varphi \right)$$

(29)

where ${\beta }_l$ is the thermal expansion coefficient and $\gamma$ is the coefficient of surface tension.

3.2.6. Definition of the Heat Source and Evaporative Energy Equation

As depicted in Figure 1, the initial vapor/liquid interface is subjected to a Gaussian distribution of the laser heat flux before deformation. Evaporation loss occurs due to the laser’s high energy density, with the laser energy calculated by Equations (30) and (31):

$$q_{Laser}=\frac{n{P}_{Laser}}{\pi R_{eff}^2}\mathit{exp}\left(\frac{-n{r}^2}{R_{eff}^2}\right)B_t$$

(30)

$$B_t=\left\{\begin{array}{c}1,t\le t_p\\ 0,t>t_p\end{array}\right.$$

(31)

where $P_{Laser}$ and $B_t$ are the laser peak power and the temporal laser distribution used to apply the effect of the pulse wave laser welding, respectively. $R_{eff}$ is the effective spot radius of the laser beam and $t_p$ is the pulse duration. The parameter n denotes the form factor for the Gaussian distribution, which is set to 2 in the present paper, as used also by Zhang et al. [34]. Moreover, in our approach, constant laser absorption is considered throughout the simulation domain. It is worth mentioning that aluminum is known to be a highly reflective material with high thermal conductivity, particularly in the spectrum commonly used in laser welding. As a result, the absorption coefficient of aluminum is relatively low compared to other materials, and the absorption of laser energy by the material remains constant within the typical operating parameters of this stationary pulsed spot laser welding. Furthermore, for the specific geometry and conditions considered in our simulation, the distribution of laser energy and absorption within the workpiece is governed by factors such as Gaussian laser beam profile, spot size, and material properties. These factors contribute to a relatively uniform distribution of absorbed laser energy within the keyhole region, supporting the assumption of constant absorption.

To apply the laser heat flux on the material surface and account for the energy loss due to evaporation, the delta function of the LS variable is utilized again. In fact, the energy loss due to evaporation is subtracted from the laser energy density, and the resulting energy is multiplied by the delta function of the LS variable to exert these effects on the vapor/liquid interface. This resulting energy is then added to the energy equation as a body force. The modified energy equation is presented in Equation (32):

$$\rho C_p\frac{\partial T}{\partial t}+\rho C_p\overrightarrow{u}.\nabla T=\nabla .\left(k\nabla T\right)+(q_{Laser}-Q_{vapor})\delta \left(\varphi \right)$$

(32)

$$Q_{vapor}=-L_V{\dot{m}}_{H-L}$$

(33)

3.3. Numerical Schemes

The entire 2D axisymmetric model was developed using COMSOL Multiphysics 5.6 within the computational fluid dynamics module to calculate the fluid flow, heat transfer, and LS transport equations. A total of 37,500 elements were used for the computational domain after conducting a mesh independency check, and the mapped typed mesh with quadrilateral elements was used with extra fine meshes as depicted in Figure 3, which were calibrated for fluid dynamic calculation. The maximum and minimum element sizes were equally set to 0.02 mm. The simulation lasted for around 17 h for 10 ms of the laser welding process, and all the calculations were done using a classical computer equipped with an Intel^® Xenon^® Gold 5118 CPU with 12 cores, 24 logical processors, and 128 GB RAM. All the time steps were chosen and set on 10 μs. A PARDISO direct solver with a preordering algorithm of Nested dissection multithreaded was used for the fluid flow calculation. The LS transport and heat transfer equations were solved using a PARDISO direct solver with the preordering algorithm of the automatic feature.

Sensitivity Analysis of the Numerical Parameters

In order to obtain reliable outcomes, the independency verification of the number of mesh elements, values of the reinitialization parameter, and the parameter controlling the interface thickness of the level-set technique were checked before the simulation process. The convergency trend and the keyhole depths were monitored as criteria for choosing the best parameters when conducting the sensitivity analysis. The analysis was done iteratively be refining simulation parameters through the steps given in

Table 3. List of the investigated tests for sensitivity analysis of the used parameters, namely, number of mesh elements (NOME), reinitialization parameter (${\mathit{\gamma }}_{\mathit{l}\mathit{s}}$), and the parameter controlling the interface thickness (${\mathit{ϵ}}_{\mathit{l}\mathit{s}}$).

Steps	Constant Values	Test Parameter	Values	Keyhole Depth
Step 1	Time step: 10 μs	NOME	16,968	4.128 mm
			24,320	4.011 mm
			37,500	3.948 mm
			48,045	3.921 mm
Step 2	${\gamma }_{ls}=3\mathrm{m}/\mathrm{s}$	${ϵ}_{ls}$	0.01 mm	Convergency error
			0.02 mm	4.058 mm
			0.03 mm	3.948 mm
			0.04 mm	3.942 mm
Step 3	${ϵ}_{ls}=0.03\mathrm{m}\mathrm{m}$	${\gamma }_{ls}$	1 m/s	4.175 mm
			3 m/s	3.948 mm
			5 m/s	3.837 mm
			7 m/s	3.839 mm

. Initially, the mesh sensitivity analysis was performed by using four different numbers of mesh elements: 16,968, 24,320, 37,500, and 48,045. The total number of mesh elements is then selected to be 37,500 since the keyhole morphology and depth showed no considerable difference with further increase in mesh elements. Then, the analysis was conducted on level-set parameters, and the most optimum values for the reinitialization parameter and the parameter controlling the interface thickness were found to be 5 m/s and 0.03 mm, respectively. The chosen parameters showed better computational efficiency and convergency trend, and further variations in the parameters did not lead to significant changes in the keyhole depth and morphology.

4. Results and Discussions

This section presents the results of the investigated cases in detail. Temperature variations within the base metal, vapor/liquid interfaces, and keyhole/molten pool shapes are examined consecutively. The first part covers the validation of the numerical model developed in this study by comparing it with a previous experimental study. Then, an overall justification is provided for the underlying physical phenomena occurring during the laser welding process, along with an exploration of the methods proposed and employed to address present interfaces and phase changes. The subsequent subsections are devoted to the impact of various laser characteristics and power wave modulations on keyhole penetration and propagation.

4.1. Accuracy Verification of Simulation Results Using Experimental Validation

The proposed model in this paper is validated by comparing it with experimental results obtained by Qin et al. [33] Specifically, the keyhole penetration depth, shape, and diameter are compared to the morphology of the keyhole observed by Qin et al. [33]. To ensure the physical comparability of the results, the simulation in this study employs the same laser characteristics as those used by Qin et al. [33]. Correspondingly, the laser power and its pulse width were set to 18 J and 3 ms. Additionally, the laser’s focal length and spot radius were both set to 60 mm and 300 µm, respectively. According to both simulation and experimental findings, the keyhole penetration depth was 3.837 mm numerically, closely aligning with the experimental value of 3.824 mm reported by Qin et al. [33]. The morphology of the keyhole propagation was also compared with experimental results, as shown in Figure 4. The spot diameter of the keyhole on the surface in the simulation and experimental result were approximately 0.937 mm and 0.936 mm, respectively. Moreover, the maximum keyhole width reported experimentally was approximately 0.407 mm. Also, the keyhole width from the numerical simulations varies between 0.310 mm and 0.51 mm. It was found that the proposed numerical approach can well predict the keyhole both in depth and diameter. Good agreements were achieved between numerical and experimental results. On the other hand, there was a slight six to twelve percent deviation between the results achieved for the keyhole width, and this is attributable to several factors. Firstly, the material properties used in the simulation might not have been the same as the material properties presented in the experimental article. Secondly, the multiple reflections of the laser beam were neglected in the simulation. Hence, laser attenuation does not occur after each reflection, which keeps the laser density at a reasonable magnitude, leading to wider keyholes. Moreover, the metallic vapor in this paper was assumed to be transparent to the laser beam, and the plasma was also neglected. Together, these factors therefore affect the input energy arriving at the material surface, which is believed to be the main reason for the slight deviation between experimental and numerical results. Nevertheless, the present model successfully correlates with the experimental results in predicting the keyhole diameter on the surface and the keyhole depth, in addition to offering acceptable magnitudes for the keyhole width.

4.2. Physical Phenomena in Laser Welding

Figure 5 and Figure 6 depict the key aspects of the physical phenomena involved, such as the keyhole, molten pool, driving forces, pressures, mushy zone, and solidus/liquidus temperature lines. A series of driving forces, including the recoil pressure, buoyancy force, hydrodynamic and hydrostatic pressures, gravity, and surface tension effects (Marangoni and curvature effect), impose fluctuations of the keyhole and its propagation throughout the laser welding process. Figure 5 depicts the driving forces acting upon the keyhole and the molten pool. In this regard, the recoil pressure, which is colored green, aims to open the keyhole, while the hydrostatic pressure (brown) and surface tension forces (purple) try to collapse and close the keyhole [38]. The buoyancy effect due to natural convection within the molten pool and the gravity are specified in the figure in red and blue, respectively. The competition between these forces is the main reason for the keyhole instabilities, specifically for regions with higher temperature gradients, such as the liquid/vapor and solid/liquid interfaces. Moreover, the solidus/liquidus temperature contour line and the mushy zone (partially-molten material) are also visible in the picture.

Figure 6 presents the laser welding process stages in greater detail, using velocity contours in which black and white arrows indicate the material’s flow direction in fluid 1 (keyhole and above the surface) and molten pool sections, respectively. Initially, the fusion begins when the laser is turned on, forming a solid/liquid interface as the temperature has not yet exceeded the vaporization threshold, which is evident within the first 30 µs. After 60 µs, the temperature rises above the vaporization point, forming a depression due to the recoil pressure following the localized vaporization. The pressure then pushes the liquid up and out of the keyhole, helping its penetration inside the material. This is corroborated by following the keyhole formation and propagation, as shown in the subsequent times after 60 µs of the process where the keyhole starts to form. As the keyhole penetrates deeper into the material, the temperature gradients, surface tension effects, and other driving forces within the molten pool and keyhole intensify, leading to increased instabilities and fluctuations in the keyhole walls. These instabilities and fluctuations are more profound around the vapor/liquid interface, where the temperature gradients are amplified. This is observed after 3000 µs, where more arrows pile up around the fluctuating interfaces.

4.3. Analyzing the Impact of Laser Characteristics on the Morphology of the Keyhole

This section presents the impacts of the laser spot radius, laser frequency, and laser power on the keyhole penetration depth, and its propagation. Details of the laser characteristics used in this section are given for cases LC1 to LC10 in Table 2. For this section of the results, a rectangular pulse shape was employed for the laser density profile, as depicted in Figure 2a.

4.3.1. Effect of Spot Radius

This section examines the impact of the spot radius on the keyhole propagation and depth by analyzing four different spot radii under a constant laser power of 6 kW. Figure 7 illustrates the keyhole morphology for each spot radius during one pulse of laser welding. The results demonstrate that the more the spot radius is enhanced, the smaller the keyhole penetration depth becomes. It is observed that larger spot radii intensify the effect of the surface tension, reducing the tendency for laser penetration in the vertical direction while increasing it towards the sides. Consequently, smaller spot radii exhibit deeper keyholes. This behavior can be attributed to the fact that smaller spot radii concentrate the laser beam on a smaller area with a higher power, resulting in a faster temperature rise at the attack surface and quicker keyhole formation. The faster formation of the keyhole corresponds to the earlier dominance of the recoil pressure, which accelerates the keyhole propagation rate.

4.3.2. Impact of Laser Frequency

This section presents the impact of the laser frequency on the laser welding process through details given in Table 2 for cases LC5 to LC7. Simulations were conducted under 6 kW laser power using three laser rectangular pulses, with frequencies of 50 Hz, 100 Hz, and 150 Hz, to investigate their impact on the keyhole penetration depth and its propagation. A spot radius of 300 µm was chosen for this analysis due to its ability to achieve a deeper keyhole and a smoother surface melt ejection than with other spot radii. The pulse width is kept constant, while an increase in the laser frequency results in shorter pulse periods. Figure 8 illustrates that the keyhole penetration depth is gradually enhanced with higher laser frequencies. Notably, even though the total laser application time is shorter at 150 Hz than at 50 Hz, the keyhole penetration depth is greater at 150 Hz. This characteristic may be advantageous for high-speed production systems. However, as depicted in Figure 8 and Figure 9, higher laser frequencies lead to heightened keyhole instability and wall fluctuations. Figure 9 shows that at a laser frequency of 150 Hz, the keyhole becomes more unstable with steeper wall slopes after two pulse periods, enhancing the risk of sudden keyhole collapse and the formation of internal porosities.

4.3.3. Impact of Laser Power

Based on the findings of the previous section for cases LC1 to LC7, a laser frequency of 100 Hz was selected as it offers a greater keyhole penetration depth than 50 Hz while exhibiting lower keyhole instabilities than 150 Hz. This section focuses on investigating the impact of the laser power on the keyhole penetration depth and instabilities under a constant frequency and pulse width during one laser welding pulse with a rectangular shape, as presented in greater detail in Table 2 for cases LC8 to LC10. Figure 10 illustrates a significant enhancement in keyhole penetration depth as the laser power is raised from 2 kW to 6 kW. This can be attributed to the accelerated temperature rise at the material surface resulting from a higher laser power, which facilitates a faster keyhole formation and propagation, and ultimately leads to deeper keyhole penetration depths.

4.4. Analyzing the Impact of Modulated Wave Welding on the Morphology of the Keyhole

In this section, the impact of MW welding on keyhole propagation is investigated through the manipulation of pulse numbers, pulse width, and pulse shapes. Details of the laser characteristics used for this section are given for cases MW1 to MW17 in Table 2. The laser density profile employed for this section of the results uses various pulse shapes, as depicted in Figure 2.

4.4.1. Impact of Pulse Width

The effect of increasing the pulse width on the keyhole morphology is examined by using four different pulse widths (0.5 ms, 1 ms, 2 ms, and 3 ms) during one pulse of the laser with a rectangular shape, as specified for cases MW1 to MW4 in Table 2. Consistent with the previous sections, a laser peak power of 6 kW and a spot radius of 300 µm were chosen. Figure 11 illustrates that as the pulse width is increased, the keyhole penetration depth also improved. In fact, with higher pulse widths, the laser is irradiated on the surface over longer durations. Hence, the temperature rise in the material increases faster, leading to a faster keyhole initiation and propagation.

4.4.2. Impact of Pulse Number

Here, the effects of different numbers (2, 6, 10, 14, and 18) of rectangular pulses over a period of 0.01 s of the laser welding process are compared. The details of the investigated cases (MW5 to MW9) for this parameter are given in Table 2. The number of pulses was increased while the total welding time and pulse width were held constant at 0.01 s and 0.5 ms, respectively. The laser power and frequency were also maintained at 4 kW and 100 Hz. As the number of pulses is enhanced, the period of each pulse is reduced, resulting in shorter intervals between laser pulses on the material surface. The keyhole morphology is then analyzed to understand the impact of the pulse number on the keyhole dynamics. Figure 12 demonstrates that increasing the pulse number leads to a greater keyhole penetration depth within the 0.01 s welding duration. This is due to the extended interaction time between the laser and the material. However, compared to single pulse cases with higher pulse widths, although the final keyhole shape at the end of each test has a wider and more curved shape, more fluctuations and instabilities are observed in the keyhole walls in the middle and initial regions of the penetration. This is due to the repeated on-and-off nature of the laser when using multiple pulses during the welding process. When the laser is turned off, the influence of the recoil pressure is significantly diminished, leading to a greater tendency for the keyhole to collapse. This is primarily driven by the effects of surface tension effect, hydrostatic and dynamic pressure, and other driving forces that contribute to closing the keyhole and causing fluctuations in the vapor/liquid interface. Hence, greater instabilities and wall fluctuations are observed when using multiple pulses through pulse number augmentation.

4.4.3. Impact of Pulse Shape

In this section, the impact of power modulation is investigated for cases MW10 to MW17, using distinct pulse shapes, as depicted in Figure 2 and Table 2. A constant average laser energy of 20 J was applied for all cases, with a welding duration of 10 ms, to ensure comparability. The results are compared with those of CW welding. Figure 13 shows the keyhole propagation morphology for different pulse shapes, including single-shape pulses (a–f) and multi-shape pulses (g–i). Among these, rectangular, trapezium (type 2 and type 1), and rectangular–triangular pulse shapes exhibit the highest keyhole penetration depths, 2.65 mm, 2.42 mm, 2.37 mm, and 2.13 mm, respectively. Rectangular–rectangular, variant–rectangular, and rectangular–trapezium pulse shapes offer intermediate depths of 1.9, 1.85, and 1.77 mm, respectively. CW welding demonstrates the lowest penetration depth of 1.51 mm, along with a smaller surface hole diameter. Notably, the variant–rectangular pulse shape produces a more stable and cylindrical keyhole, with a smaller keyhole depth/width ratio than other pulse shapes.

4.5. Temperature Variations within the Base Metal

This section focuses on measuring and comparing the maximum temperatures within the base metal to gain insights into temperature variations in the laser welding process.

Figure 14, Figure 15 and Figure 16 depict the maximum temperature variations within the different base metal cases. Figure 14a–c displays the temperature variations, considering distinct spot radii, laser frequencies, and laser power, respectively. Figure 14a shows that the incidence of the laser beam on the material surface causes a sudden temperature increase, triggering a fusion phase change as the temperature exceeds the melting temperature. Subsequently, the temperature reaches approximately 2743 K, which is the vaporization threshold, and vaporization becomes an additional phase change added to the process alongside fusion. During a single pulse, the temperature changes smoothly upon reaching the vaporization point, reflecting dominant latent heat transfer and leading to a nearly constant temperature throughout the phase transition. After 2 ms (1 pulse), the laser shuts down, resulting in a drop in the maximum temperature within the base metal. The material begins to cool down until it reaches the solidification temperature, indicated by a dashed circle, leading to the solidification phase change. The temperature then decreases smoothly throughout the solidification process since the material loses its heat in a latent heat form. Further away from the solidification temperature interval, the material experiences steeper temperature drops as heat is lost in a sensible form until it reaches a stable temperature almost equal to ambient. Figure 14a also shows that the temperature reaches its vaporization point faster for smaller spot radii (LC1). This leads to faster keyhole formation, resulting in faster and deeper penetration depths, as corroborated in Figure 7. In Figure 14b, the impact of different laser frequencies on temperature variations is examined using three pulses with reduced pulse periods, revealing similar trends for temperature variation as observed in Figure 14a. For a laser frequency of 50 Hz, the maximum temperature diminishes considerably after one pulse, reaching temperatures close to ambient. In contrast, frequencies of 150 Hz and 100 Hz exhibit higher temperatures at the end of the first pulses. This difference is attributable to the longer pulse periods in these cases, which allow the material to cool down more between each pulse. Another noteworthy feature, indicated by different line arrows in Figure 14b, is the gradual increase in maximum temperature within the domain after each consecutive pulse. This is attributed to the overall temperature escalation in the base metal as the number of pulses increases. The phenomenon is more pronounced at 150 Hz due to its shorter pulse period, which limits the cooling of the material between pulses. Figure 14c shows that higher laser powers accelerate keyhole formation and propagation by allowing to reach the vaporization threshold faster. Moreover, the solidification process is prolonged at higher powers due to increased fusion throughout the material.

Figure 15a,b depict the maximum temperature change for different pulse widths and pulse numbers. Figure 15a shows the maximum temperature variation, while the pulse width enhanced from 0.5 ms to 3 ms. It is seen that longer solidification periods are observed as the pulse width increases. This might lead to higher keyhole fluctuations due to the impact of solidification on the interfaces, intensifying the surface tension effect. Moreover, with an extended pulse width, the vaporization process is prolonged, resulting in deeper keyholes, as corroborated in Figure 11. Another visible feature from Figure 15a is that for greater pulse widths, the final temperature of the material at the end of one pulse period is higher. This is due to the overall intensification in the material temperature thanks to more extended laser/material interactions induced by longer pulse widths. Figure 15b illustrates the impact of an increase in the number of pulses through a fixed welding duration of 0.01 s. As can be seen, the more the pulse number is raised, the more the period of pulses is reduced. Hence, the material loses its chance to cool down significantly between pulses. Therefore, as the number of pulses enhances, the maximum temperature at the end of each pulse period is increased, up to the point where no reasonable temperature loss is seen between pulses, as shown for case MW9, which applies 18 pulses. The temperature variations are then similar to what is observed in continuous wave laser welding.

Figure 16a,b depict the impact of various pulse shapes on the variations in maximum temperature within the base metal. Figure 16a reveals that for continuous wave welding and the variant–rectangular pulse shape with ascending rectangular steps, the maximum temperature variations are approximately the same after the vaporization point is reached. There are descending temperature variations for MW10, MW13, MW11, and MW12, which correspond to rectangular, trapezium (type two/one), and triangular pulse shapes, respectively. This is due to a gradual power mitigation from 4 kW to 0 at the end of their pulse widths, 5 ms, 8 ms, 8 ms, and 10 ms, respectively. The temperature drop when using a rectangular pulse shape (MW10) is faster at the end of one pulse as it experiences a sudden power drop from 4 kW to 0. On the other hand, for trapezium (type 1: MW11 and type 2: MW13) and triangular pulse shapes (MW12), the temperature drops with smaller slopes since the laser power is decreased less rapidly but rather gradually than what is observed for MW10. There are no significant differences in the maximum temperature variations within the material when using multi-shape pulses, as can be seen in Figure 16b.

5. Conclusions

This paper presents a novel numerical approach utilizing COMSOL Multiphysics to simulate pulsed and continuous laser welding of aluminum in a 2D axisymmetric configuration. The study investigates the effects of laser parameters such as the spot radius, laser power, and laser frequency on keyhole penetration depth, propagation, and instabilities. Moreover, the impact of power wave modulation is also examined on keyhole dynamics through variations in the pulse width, pulse number, and pulse shape. The results are compared with continuous welding, and the following findings are observed.

• The more the spot radius is enhanced, the smaller the keyhole penetration depth, and the more intense the melt ejection. A reduction of over 80% in the keyhole penetration depth is observed with an increase in the spot radius.
• As the laser frequency increases, the keyhole wall instabilities and the tendency of the keyhole to collapse are amplified while the keyhole penetration depth is increased to some extent.
• With an increase in laser power from 2 kW to 6 kW, the keyhole penetration depth is improved by more than 80%.
• Extending the pulse width from 0.5 ms to 3 ms leads to an increase of over 80% in the keyhole penetration depth. Moreover, the keyhole wall becomes more unstable as pulse width is extended.
• If the welding duration is maintained at 0.01 s, the keyhole penetration depth increases significantly when using higher pulse numbers. However, more keyhole fluctuations and instabilities are observed due to multiple laser on-and-offs.
• The rectangular pulse shape has the greatest keyhole penetration depth among various pulse shapes, while variant–rectangular pulse shapes and triangular pulse shapes produce more keyhole stability with smaller depth/width ratios.
• At the end of the welding process, higher temperatures within the base metal, achieved during CW laser welding, do not necessarily correspond to deeper keyholes and welding efficiency.