Numerical Analysis of the Partial Penetration High Power Laser Beam Welding of Thick Sheets at High Process Speeds

: The present work is devoted to the numerical analysis of the high-power laser beam welding of thick sheets at different welding speeds. A three-dimensional transient multi-physics numerical model is developed, allowing for the prediction of the keyhole geometry and the ﬁnal penetration depth. Two ray tracing algorithms are implemented and compared, namely a standard ray tracing approach and an approach using a virtual mesh reﬁnement for a more accurate calculation of the reﬂection point. Both algorithms are found to provide sufﬁcient accuracy for the prediction of the keyhole depth during laser beam welding with process speeds of up to 1.5mmin − 1 . However, with the standard algorithm, the penetration depth is underestimated by the model for a process speed of 2.5mmin − 1 due to a trapping effect of the laser energy in the top region. In contrast, the virtually reﬁned ray tracing approach results in high accuracy results for process speeds of both 1.5mmin − 1 and 2.5mmin − 1 . A detailed study on the trapping effect is provided, accompanied by a benchmark including a predeﬁned keyhole geometry with typical characteristics for the high-power laser beam welding of thick plates at high process speed, such as deep keyhole, inclined front keyhole wall, and a hump.


Introduction
The high degree of focusability and radiation intensity of the laser beam lead to the formation of a cavity within the molten pool, called a keyhole. A cavity-enhanced optical absorption, known nowadays in laser beam welding (LBW) as the keyhole mode technique, was recognized for its high energy efficiency about 50 years ago [1]. Within the keyhole, multiple-reflection and correspondingly multiple-absorption events are observed, resulting in the enhanced optical absorption. Hereby, the total amount of absorbed energy, determined mainly by the keyhole geometry, increases, allowing for the joining of highthickness components. More recently, the single-pass LBW of thick weld samples, with reachable thickness of up to 50 mm [2], was enabled by the keyhole mode technique and the rapid development of modern laser systems with available laser power of up to 100 kW [3][4][5]. Furthermore, the unique and well-known technical advantages offered by the LBW process, e.g., low distortion and narrow heat affected zone (HAZ), allowed for a drastic increase in the industrial productivity. Hence, at the present time, the LBW process finds application in the shipbuilding and aerospace industries, in the production of thick-walled structures, such as pipelines for the oil and gas industry, or even in the manufacturing of vacuum vessels and correction coils for the International Thermonuclear Experimental Reactor (ITER) [6][7][8].
Despite the advantages mentioned above, aiming at joining specimens of high thickness by the LBW process, especially at higher process speeds, remains challenging due to the formation of welding defects [9], such as porosity [10], spatter formation [11,12], sagging of liquid metal [13,14], hot cracking [15][16][17], etc. As suggested by earlier and recent research results, the keyhole and molten pool characteristics are decisive for the final weld seam quality. In other words, a stable keyhole and molten pool result in fewer welding defects and enhanced mechanical properties of the welded component.
In order to gain more understanding about the keyhole and molten pool dynamics and more precisely how these correlate with the formed defects, several experimental attempts have been published in the literature. Monitoring and detecting systems, based on different measuring techniques, e.g., acoustic [18], infrared and imaging signals [19], high-speed imaging [20], etc., have been utilized for the collection of process data. By making use of the collected data, new insights into the process can be obtained. For example, a three-dimensional weld pool surface reconstruction has been achieved using a high-speed camera and a dot matrix pattern laser [21,22]. The keyhole size and the molten pool edges has been estimated with a coaxial monitoring system in [23,24]. A porosity prediction in real-time experimentation based on a visual monitoring system has been proposed in [25]. Due to recent advancements, a more sophisticated measuring and visualization method using high-speed synchrotron X-ray imaging to monitor the keyhole shape inside the molten pool has been developed by the Joining and Welding Research Institute of the Osaka University [26,27]. In recent publications, this method has been combined with a simultaneous monitoring of the absolute energy absorption by omnidirectional backscattered laser intensity. This enables the real-time visualization of the relationship between the laser absorption and the keyhole depth [28].
Although the experimental methods mentioned above can be used to ascertain the keyhole and molten pool dynamics, their application remains limited by the highly expensive and bulky equipment, as well as by the relatively low image resolution. Strictly speaking, no insightful and detailed information, such as transient flow pattern, velocity or temperature distributions, can be gained, making the study of the LBW process very difficult.
Due to the advancement of computational technologies in recent decades, numerical modeling has become an established research tool, allowing for the estimation of important process characteristics, such as molten pool shape and thermal cycles. Early numerical models of the LBW process focuses entirely on the conduction mode welding, thus neglecting important physical aspects of the beam-matter interaction, such as multiple reflection, evaporation, and free surface deformation. Nonetheless, it is worth noting that depending on the aim of the study, even simple heat conduction models can be sufficient, e.g., for the prediction of the fusion zone with a concentrated and uniformly distributed heat source [29][30][31] or with a non-concentrated, non-uniform heat source distribution [32,33]. On the other hand, the study of complex welding defects, such as pore formation and hot cracking, demands a detailed description of the underlying physics, especially of the beam-matter interaction and the thermo-fluid dynamics. One of the very first attempts to estimate the two-dimensional keyhole shape and the corresponding energy distribution on the keyhole wall was made by considering multiple reflections in a conical keyhole with an averaged inclination [34]. Since then, several advanced numerical models calculating the multiple reflections using a so-called ray tracing technique have been developed by the welding community. The ray tracing methodology enabled the precise prediction of the energy distribution within the keyhole by calculating the current location of each sub-ray and the location as well as the direction of its reflection. Basically, the advanced models can be subdivided into two groups according to the number of phases considered in them, namely into two-phase and three-phase models. Hereby, the three-phase models take into account the solid, liquid and gas phases, e.g., References [35][36][37][38][39] and the two-phase model consider only the solid and the liquid phase, e.g., References [40][41][42]. Note that some of the two-phase models include the vapor-induced effects empirically, e.g., References [43,44].
An overview of the available two-phase and three-phase models with a detailed description of the considered physical phenomena can be found in [45].
A critical review of the available LBW numerical models shows that a reliable model allowing for the study of the keyhole and molten pool dynamics must include an accurate description of the temporal and spatial energy distribution on the keyhole wall, e.g., by utilizing a ray tracing algorithm. Furthermore, the two-phase models considering the vapor-induced effects empirically seem to offer the best balance between computational intensity and realistic results. To the best of the authors' knowledge, the majority of the available LBW models concentrate either on the study of thin sheets welding with a thickness below 5 mm, at high process speeds, or on thick sheets welding at low processing speeds, typically less than 1.5 m min −1 . However, the great potential for more effective joining of high thickness components with the LBW process seems to be impeded by the occurrence of untypical defect formation, especially at higher process speeds. Hence, the need for a practical and at the same time reliable numerical model for the LBW of thick sheets at high process speeds arises.
The main objective of the present study is to develop a suitable approach, allowing for the accurate prediction of the keyhole and molten pool dynamics during LBW of thick sheets at high process speeds. The focus of the study is on the accuracy of the chosen ray tracing approach. This is directly related to the energy distribution on the keyhole wall, thus having the strongest impact on the keyhole and molten pool dynamics. A three-dimensional transient multi-physics numerical model is developed and tested. Two ray tracing algorithms are implemented, namely a standard ray tracing approach and an approach using a virtual mesh refinement for more accurate calculation of the reflection point. The results obtained with both ray tracing algorithms, such as penetration depth and process efficiency, are compared with experimental measurements. Finally, best practices are obtained based on the limitations of the currently available LBW models with regard to their application to the study of welding thick sheets at high process speeds.

Materials
Unalloyed steel sheets S355J2+N of 12 mm thickness were utilized in the welding experiments. The dimensions of the sheets were 175 mm × 100 mm × 12 mm. The corresponding chemical composition was measured with spectral analysis and is given in Table 1.

Experiments
All welds produced in the experiments were bead-on-plate welds performed with a 16 kW disc laser Trumpf 16002. A schema of the experimental setup is shown in Figure 1; the process parameters are summarized in Table 2. Macro sections have been extracted from the middle region of the weld seam marked in Figure 1. From these, several metallographic cross-sections have been prepared with a 2% nital etching, which subsequently were used for the validation of the numerically obtained results.

Numerical Modeling
A three-dimensional thermo-fluid dynamics model tracking the free surface deformation by the volume of fluid (VOF) technique was developed for the study of the high power LBW process at high process speeds. The multi-physics model is based on several previous works with some further improvements and adaptions. More details can be obtained from the authors' previous works [46][47][48][49][50]. Note that only the main physical features of the welding process are formulated concisely in the manuscript; thus the emphasis is on the improvements in the model, especially on the ray tracing approach.

Assumptions
In recent years, the numerical modeling of many industrial and scientific problems has become an inseparable part of their research due to the rapid increase in computational capacity. However, the precise mathematical description of the complex physics behind the high-power LBW process, including numerous strongly coupled, highly-nonlinear interactions between the laser radiation, the vapor phase, the molten metal, and the solid material, remains challenging. Thus, several assumptions and simplifications need to be made in the simulation procedure to allow for feasible computational times. The main assumptions made in the model are the following: • The molten metal and the gas phases are assumed to be Newtonian and incompressible fluids.

•
The flow regimes in the model are considered to be laminar.
• The Boussinesq approximation is used to model the impact of the density deviation on the flow.

Governing Equations
The governing equations describing the multi-physics model in a fixed Cartesian coordinate system are summarized below.
The VOF technique is applied to track the transient deformation of the molten pool free surface and solidified weld seam profile.
• Volume fraction conservation where α vol steel denotes the volume fraction of the steel phase in a control volume and is the fluid velocity vector [51].
• Mass conservation where ρ mix is the volume-fraction-averaged density, t is the time, p is the fluid pressure, µ mix is the volume-fraction-averaged dynamic viscosity, #» g is the gravitational acceleration vector, and #» S m is the source term. The source term takes into account the thermal buoyancy due to the variations of the density of the steel with temperature [52]; the deceleration of the flow in the mushy zone, which is related to the inverse of the size of the interdendritic structure [53,54]; the effects of surface tension along the steel-air interface [55]; the evaporation-induced recoil pressure [56]; and the vapor-induced effects, such as stagnation pressure and shear stress on the keyhole surface [43].
• Energy conservation where H mix is the enthalpy, λ mix is the heat conductivity, and S e is the source term. The source term takes into account the laser heat flux density, the convective and radiative heat transfer, the evaporation loss and the recondensation. Note that the outward as well as the inward convective and radiative heat fluxes due to the high-temperature metal vapor, reaching temperature of up to 6000 K, are considered according to [44]. In the present work the range of action for the vapor-induced secondary heat effects was adapted according to the geometrical dimensions of preliminary obtained metallographic cross-sections.

Raytracing Algorithms
A precise description of the laser energy distribution on the keyhole wall is crucial for an accurate calculation of the transient keyhole and weld pool dynamics. A physics-based self-consistent ray-tracing algorithm, rather than an empirical volumetric heat source, is commonly preferred for the description of the multiple reflections and the Fresnel absorptions considered as the key physical aspects of the beam-matter interaction in the LBW process. In the authors' previous studies, a ray-tracing algorithm, based on the work presented in [40], was implemented in the commercial finite volume method software Ansys Fluent via user-defined functions. Instead of calculating the radiative transfer Equation [57] or the path of Lagrangian photonic parcels [58], the laser beam is divided into 755 sub-rays, which results from discretizing the laser spot in the focal plane by 31 × 31 sub-regions, and each sub-ray is given its own location-dependent energy density and initial incidence angle according to the laser beam profile. Hereby, the laser heat flux density is assumed to have a Gaussian-like axissymetric distribution according to [43]. The temporal and spatial distribution of the laser energy on the keyhole wall can be obtained after calculating the multiple reflections and the Fresnel absorptions of all sub-rays [59].
Since the keyhole profile calculated with the VOF method is not explicit, it may lead to difficulties in determining the exact reflection position of the sub-ray geometrically. Therefore, a compromised criteria is applied in the present study, identifying a cell as a "reflection cell" when the following criteria is satisfied [40]: where d ray is the distance between the cell center and the incident ray; ∆ cell is the cell size. It can be seen from Equation (5) that the accuracy of the ray-tracing algorithm is directly determined by the cell size. However, the recommended optimal cell size of 0.2 mm [60], which is comparable with the laser spot radius of 0.25 mm, has to be chosen for an affordable computational time. This leads to certain inaccuracies, especially when the keyhole front wall is nearly parallel to the laser beam and fluctuates with an amplitude comparable to the cell size, as shown in Figure 2a, where the gray regions represent the steel phase within the cells. The deviations may become unacceptable for high-power LBW of thick sheets at high process speeds due to the higher energy densities, leading to more severe fluctuation of the keyhole wall. Therefore, a local virtual grid refinement algorithm was developed to improve the accuracy of the ray-tracing algorithm [61] without a noticeable increase in the computational intensity; see Figure 2b. By assuming the free surface to be a piecewise linear interface, the distance between the cell center and the free surface within the cell, D f ree (see Figure 3) can be determined uniquely by using the volume fraction of the steel phase α vol steel and its gradient, #» n , given as #» n = ∇α vol steel |α vol steel | = n x , n y , n z [62]. Taking the cell center as the origin and substitut-ing n 1 , n 2 , n 3 as the maximum, median, and minimum values from n x , n y , n z , five scalar conditions are obtained, as shown in Figure 3. These can be expressed as: D f ree can be calculated analytically or numerically by using Equation (6). Hereby, the free surface within the cell is approximated as: n x x + n y y + n z z = D f ree , and once D f ree is known, it can be used to determine the virtual cells lying on the free surface.
In the first step of the ray-tracing algorithm, all potential reflection cells are selected by applying the condition defined in Equation (5) and using ∆ cell = 0.2 mm. Subsequently, the selected cells are further divided into virtual cells with a typical cell size of 0.05 mm. The virtual cells on the free surface are identified by Equation (7) and marked in yellow as seen Figure 2b. Note that the cells in Figure 2b marked in gray are the cells excluded from the calculation of the reflection point. In the last step, the reflection point is determined among the selected cells by repeating step one with the virtual mesh cell size ∆ cell virtual = 0.05 mm.

Boundary Conditions
Following the basic principles of fluid dynamics and assuming the air phase to be inviscid (µ air = 0) and incompressible, two scalar conditions for the pressure and viscous stress on the steel-air interface are defined: Note that the surface unit normal is directed into the interior of the steel phase. The energy boundary condition on the steel-air interface, considering the multiple Fresnel absorption, heat convection, thermal radiation, evaporation, and recondensation, is expressed as: The boundaries of the air-phase domain were set as pressure outlets and on the bottom of the steel-phase domain, convective heat transfer was considered. The simulation domain was smaller than the real steel sheets to optimize the computational time. Hence, a proper treatment of the boundaries was ensured by utilizing the continuity boundary condition [63]. The boundary conditions are also shown in Figure 4.

Material Properties
Temperature-dependent material properties were implemented for the simulation of the LBW process. The maximum reachable temperature in the simulation was set to 3400 K. The base material was modeled as ferritic phase, whereby for the temperatures above the austenitization temperature, the properties of the austenite phase were used. The phase-specific properties were taken from the literature [64][65][66]; if they were not available, the values for pure iron were taken, due to the close chemical composition (Fe∼98%). The solid-solid phase transformation (ferrite-austenite) was included to the specific heat and the austenite-martensite phase transformation was neglected due to its relatively small amount of latent heat. The thermo-physical material properties are shown in Figure 5. An averaged value of the steel density was calculated for the temperature range of interest between 1200 K and 2800 K, giving ρ steel = 7060 kg m −3 .

Numerical Setup
The computational domain used in the present study had dimensions of 18 mm in length, 8 mm in width, and 12 mm in thickness; see Figure 4. The domain is uniformly meshed by hexahedral cells of 0.2 mm, resulting in a total amount of 189,000 control volumes. An air layer between 0 mm ≤ z ≤ 2.0 mm was defined above the steel sheet (2.0 mm ≤ z ≤ 14.0 mm), allowing the tracking of the steel-air interface by the VOF method.
All governing equations were solved with the commercial finite volume method software ANSYS Fluent. A second order upwind scheme was used for the spatial discretization of the momentum and energy conservation equations, and a first-order implicit formulation was applied for the discretization of the transient terms. The pressure-velocity coupling was realized by the PISO scheme and the steel-air interface was reconstructed by the Geo-Reconstruct method.
A high-performance computing cluster with 88 CPU cores at the Bundesanstalt für Materialforschung und -prüfung (BAM) was used for the computation. The averaged computing time was approximately 24 h for 0.25 s real process time.

Process Efficiency and Drilling Speeds
First, the process efficiencies for both welding parameters, given in Table 2, with the different ray tracing algorithms described above are calculated and compared. However, although different welding parameters were used, the line energy was kept equal. A closer look at Figure 6 shows that there are much less fluctuations in the amount of absorbed energy predicted by the standard ray tracing method compared to the virtually improved algorithm, for both low and high process speeds. There were almost constant maximum values of the absorbed laser energy during the simulation with the standard algorithm; see Figure 6a,c. In contrast, the simulation with the virtual mesh refinement does not have clearly pronounced maximum value for the absorbed energy; see Figure 6b,d. Strictly speaking, this means that in the calculations with the virtually refined ray tracing method, more sub-rays are either reaching the keyhole bottom or have been reflected out of it, leading to fluctuations around the averaged value of the absorbed laser energy. Hence, the averaged values of the absorbed amount of laser energy obtained with the improved ray tracing method are slightly lower than those predicted by the standard algorithm. The averaged values at the lower process speed of 1.5 m min −1 are 83% and 76% for the standard and virtually refined method, respectively. At the higher process speed of 2.5 m min −1 , the differences are even smaller, with 83% for the standard and 81% for the virtually refined approach. Note that all calculated process efficiencies lie within the experimental range of measurements at similar process conditions; see, e.g., [67]. Thus, no conclusive statement regarding the physicality and accuracy of the ray tracing algorithms can be given. In the next step, the drilling curves and speeds predicted by the standard and virtually refined ray tracing algorithms are compared and analysed. As seen in Figure 7a, at lower process speed, both algorithms deliver similar results. The penetration depth reached is about 7.7 mm and 8.5 mm for the standard and virtually refined methods, respectively. Furthermore, the drilling speed in the first 50 ms of the simulation, estimated from the slopes of the drilling curves, seems to be similar for both ray tracing methods (94 m s −1 vs. 124 m s −1 ). On the other hand, a closer look at the results obtained for LBW at higher process speeds (Figure 7b) shows much bigger differences in the final penetration depth as well as in the drilling speed. The deviation of the penetration depth is about 3 mm, and the drilling speed is almost two times higher in the simulation utilizing the virtually refined technique (114 m s −1 vs. 196 m s −1 ). This result is surprising, since both welding parameter sets had the same line energies, and the differences of the calculated process efficiencies were negligible from a practical point of view. In the last step, the accuracy of the numerically obtained penetration depths from both algorithms is validated by comparing them to experimental measurements; see Figure 8. It can be seen from Figure 8a that in the case of modeling LBW at comparatively lower process speed of 1.5 m min −1 , both ray tracing algorithms deliver results that are close to those from the experiments. It is worth mentioning that although the penetration depth obtained with the standard method is lower than the value predicted by the improved ray tracing algorithm, both values lie within the tolerance range of the experiments. Hereby, the tolerances of the simulation are determined by the half size of the control volumes of 0.1 mm. The experimental tolerances on the other hand are obtained by comparing three metallographic cross-sections extracted from the quasi-steady state region of the weld seam, according to Figure 1. In Figure 8b,c exemplary metallographic cross-sections of 12 mm thick unalloyed steel sheets laser beam welded at 1.5 m min −1 and 2.5 m min −1 , respectively, are shown. As seen in Figure 8a, at the higher process speed of 2.5 m min −1 , significant differences emerge between the numerically predicted by the standard method and experimentally obtained penetration depths. However, the results obtained with the virtually refined ray tracing algorithm show very good agreement with the measured penetration depth. Hence, it can be concluded that the standard ray tracing technique is well-suited for the computation of the LBW process of thick sheets at lower process speeds. However, the virtually refined ray tracing algorithm seems to provide sufficient accuracy at both low and high welding speed.

Correlation between the Energy Distributions in the Keyhole and the Penetration Depth
In order to reveal the physical reasons for the different accuracy of the ray tracing algorithms regarding the predicted penetration depths, the energy distribution in the keyhole is studied in detail for both methods. In Figure 9, the location of the maximum laser energy during LBW at 2.5 m min −1 is exemplarily plotted in the time interval between 0.25 s and 0.26 s for the standard and virtually refined algorithm. Note that the maximum value of the absorbed laser energy is used for the analysis since it provides comprehensive information about the keyhole geometry as well. For example, according to the basic principles of the LBW process, the ablation of the front keyhole wall is directly connected with the movement of the maximum laser energy location from the top to the bottom region of the keyhole. Once the the energy is concentrated on the keyhole bottom, the keyhole depth increases, allowing for the deep-penetration LBW [68]. Thus, the ablation process of the front keyhole wall is crucial for the LBW process, especially during welding of thick plates, when it is most pronounced.  The in-depth analysis of the graphs seen in Figure 9 shows that the maximum laser energy location follows the typical pattern for the LBW process, by first ablating the front keyhole wall and subsequently concentrating on the keyhole bottom, for both ray tracing methods. For instance, during the simulation with the standard algorithm, the ablation process is easily recognized in the time intervals between 0.25 s and 0.253 s, as well as between 0.256 s and 0.26 s; see Figure 9a. However, the standard algorithm seems to frequently trap the maximum laser energy in the top region, e.g., in the interval between 0.253 s and 0.256 s. Logically, due to this trapping effect, the ablation time of the front wall increases, leading directly to a decrease in the duration time of the maximum energy on the keyhole bottom, thus resulting in a lower penetration depth. In contrast, the trapping effect is not apparent when using the virtually refined method, providing greater accuracy. Additionally, to better visualize the physicality of the ray tracing algorithm, a benchmark was developed and utilized. The benchmark includes a predefined keyhole geometry with typical characteristics for the high-power LBW of thick plates at high process speed, such as deep and narrow profile, an inclined front keyhole wall, and a hump thereon, i.e., a dynamic fluctuation. Both ray tracing algorithms are applied to the predefined keyhole geometry, as shown in Figure 10. With a closer look at Figure 10a, it can be easily seen that the standard algorithm traps most of the laser energy in the hump region. Thus, only a small part of the total laser energy reaches the bottom, resulting in the insufficient penetration depth. However, as seen in Figure 10b, part of the laser beam can pass the hump and directly reach the bottom of the keyhole when the virtually refined raytracing method is applied. Additionally, without the trapping effect part of the laser energy is first reflected to the rear keyhole wall and subsequently reflected to the keyhole bottom. Hence, much more energy is accumulated at the bottom of the keyhole, leading to realistic predictions of the penetration depth. It is, however, worth mentioning that the increase in accuracy of the calculated energy distribution within the keyhole can also be achieved by directly using finer mesh, i.e., of the size of the virtual mesh of about 0.05 mm. However, such discretization will increase the number of control volumes by approximately 64 times when meshing the computational domain uniformly with hexahedral elements. Additionally, to ensure numerical stability, smaller time steps should be used due to the smaller cell size. Hence, the use of finer discretization remains challenging and unrealistic, even at present time. In contrast, the virtually refined ray tracing technique allows approaching the problem successively and obtaining an accurate solution by almost no increase in the computational intensity.

Conclusions
In the current work the plausibility of the well-known and widely used ray tracing algorithms is studied for the high power LBW of thick sheets at different process speeds. The main conclusions drawn based on the obtained results are summarized as follow: • A three-dimensional transient multi-physics numerical model is developed, allowing for the prediction of the keyhole geometry and the final penetration depth.
• Two ray tracing algorithms are implemented, namely a standard ray-tracing approach and an approach using a virtual mesh refinement for more accurate calculation of the reflection points. • Both algorithms provide sufficient accuracy for the prediction of the keyhole depth during LBW with process speeds of up to 1.5 m min −1 . • The standard algorithm is found to underestimate the penetration depth at process speeds of 2.5 m min −1 due to a trapping effect of the laser energy in the top region. • The virtually refined ray tracing approach results in high accuracy results for both process speeds of 1.5 m min −1 and 2.5 m min −1 .