3D Thermal Simulation of a Laser Drilling Process with Meshfree Methods

Abstract

Numerical simulation of laser drilling is rapidly gaining interest in academia and industry since this process remains one of the most important and widely-used technologies in modern manufacturing. Meshfree methods such as Smoothed Particle Hydrodynamics (SPH) have proven to be successful as a numerical tool for the computation of the heat transfer and material removal associated with a laser drilling problem. Nonetheless, the vast majority of recent developments incorporate an inconsistent SPH kernel into their thermal simulations. In this paper, several enhanced schemes are implemented to address this problem by solving the heat transfer more accurately. These meshfree schemes can provide a second-order accurate discretization of the Laplace operator and abolish the inconsistency issue of the standard SPH kernels. An efficient approach is additionally suggested to handle the associated boundary conditions, which relies on the idea of the color function and particle label. The implementation is initially validated by a 3D benchmark study and then applied for the first time to a laser drilling problem.

1. Introduction

Laser drilling is a manufacturing process where pulsing laser beams are focused on a workpiece to create popped holes (also known as thru-holes). Due to their unique efficiency and accuracy, these technologies have become the center of attraction in various fields of application such as the gas turbine, aerospace, and printed circuit board industries. The analysis of a laser drilling problem is typically governed by an energy balance, where the thermal energy provided by a focused laser beam is conducted into the workpiece and consumed to melt or vaporize the material. This description indicates that solving the heat equation for solids is of crucial concern for the analysis of laser drilling problems since thermal effects are usually the prevalent mechanism during such processes. In brief, there exists a wide range of physical phenomena (including the thermal issues), process variables, and material parameters in laser manufacturing, recognized as the main complexities for its theoretical and experimental developments.

Numerical modeling, on the other hand, is an alternative approach that can leverage the improvement potentials of laser drilling processes. Reliable computational models are valuable tools for the characterization of quantities that are impossible to predict by theory and/or extremely hard/expensive to measure in experiments. From this point of view, laser drilling problems seem to be a prime candidate for numerical simulation since they are typically involved in ultra-fast pulses and abrupt phase changes [1,2]. According to hundreds of scientific citations, the Finite Element Method (FEM) may be regarded currently as the tool-of-choice in the domain of laser processing simulations. Examples include the FEM model designed for high energy lasers by [1], the investigation of temperature fields in sheet metals [3], the modeling of the laser welding process in butt-joint specimens [4], and the simulation of fundamental laser-micromachining processes [5]. While the use of the Finite Difference Method (FDM) [2] and the Finite Volume Method (FVM) [6] can also be found in the literature, it is hard to overstate the superiority of FEM over other mesh-based manufacturing simulations.

Another class of spatial discretization is based on particles (instead of a mesh), which is found to be a particularly promising approach for applications with complex free-surfaces, phase changes, and multiple material interfaces. The simulation of laser-based manufacturing technologies in general, and laser drilling in particular, encounters all these difficulties to a considerable extent. Therefore, several researchers have attempted to investigate the performance of meshfree particle methods in laser-based manufacturing simulations. In 2008, Gross [7] introduced a first version of the Smoothed Particle Hydrodynamics (SPH) model for modeling a simple laser cutting process. Later, an SPH thermal model of the direct laser interference patterning was developed by [8]. In the article published by Kim [9], Isoparametric Finite Point Interpolation (IPIM) techniques in both weak and strong forms were applied to simulate the material removal process in laser drilling. These early studies demonstrated an enormous potential of meshfree methods in simulating laser manufacturing processes.

More recent works include the 3D laser cutting simulation by [10] using a commercial SPH package and an improved SPH formulation for heat transfer in the laser ablation of aluminum [11]. Other interesting efforts, more closely falling into the scope of this paper, are the application of the Radial Point Interpolation Method (RPIM) [12] and the Meshless Local Petrov–Galerkin (MLPG) collocation method [13] for the heat conduction in laser drilling. In [14], the authors reproduced the results of [12] by using yet another particle-based scheme. It is worth mentioning that the number of SPH developments in other fields of application is overwhelming, and numerous efforts have been made to model a variety of structural [15], CFD [16,17], and manufacturing simulations [18,19,20,21], to name a few. While a growing share of interest in meshfree simulations of laser processing can be noticed upon this literature survey, a comprehensive study of higher order SPH formulations for (3D) thermal modeling in laser-based manufacturing processes seems to be missing from the current state-of-the-art.

Thus, the ultimate goal is to develop a highly efficient and robust meshfree solver that can provide a design tool for laser-based manufacturing models including all thermal-mechanical-material coupling effects. In the pursuit of a framework with such capabilities, this work sets its sights on the development of a robust 3D thermal solver for laser drilling applications. The first application of the methods listed in

Table 1. Summary of the methods implemented in this work.

Scheme	Acronym	Order of Consistency	Year	Publication
BrookSHaw’s approximation	BSH	–	1985	[22]
Particle Strength Exchange	PSE	–	1986	[23]
Reproducing Kernel Particle Method	RKPM	1st	1995	[24]
Corrective Smoothed Particle Method	CSPM	1st	1999	[25]
Fatehi’s Meshfree Scheme	FMFS	2nd	2011	[26]
Improved CSPM	ICSPM	2nd	2017	[27]

to the heat transfer simulation of a laser drilling model is presented. Furthermore, an efficient strategy for handling the boundary issues is suggested, which relies on the idea of the color function and particle label. We also highlight the use of more advanced SPH formulations to recover first- and second-order consistency for the Laplace operator. This contribution allows for a more accurate solution of the heat equation, especially in the presence of complex boundary conditions. After validating the computer implementation in a 3D benchmark test, the proposed framework is subsequently utilized to simulate a laser drilling process.

The remainder of this paper is organized as follows. In Section 2, the relevant physics is described in some detail, including the more important of the simplifications employed by multiple research groups active in the field. Section 3 is the core of this article, where the selected meshfree methods are presented in such a way that they can be implemented in computer code with minimum effort. In Section 4, the numerical simulation framework is first validated by performing a convergence study and then applied to the laser drilling process. Concluding remarks are presented in Section 5.

2. Laser Drilling Model

A graphical illustration of the laser drilling process is displayed in Figure 1, where a laser beam is acting on a metallic material. As can be viewed in this figure and according to [28], laser drilling of cylindrical holes generally occurs through the melting and vaporization of the work material, plus the auxiliary action of the assist gas. In such a process, the focused laser beam serves as the external source providing the energy required to remove material. On the other hand, the assist gas facilitates the drilling process via an exothermic reaction (usually with oxygen) and/or evacuates molten material from the popped hole.

2.1. Problem Statement

A complete physical model of such processes is far from simple and would warrant original research as it needs to resolve multiple challenging issues including new surface generation, different time and length scales, and multiple material-phase interfaces with extreme and complex movement. In this paper, the diverse range of physical phenomena governing laser drilling processes is reduced to a thermal simulation in solids such that the present numerical solver can be applied effectively. In what follows, the details of the underlying physical model are described.

In our simplified model, a static laser beam source with a $\mathsf{\mu }$s to ms pulse duration is considered, while the assist gas is not modeled in the thermal analysis. The single-pulsed laser used here has a pulse duration of 150 $\mathsf{\mu }$s and a peak power of 100 W exerted on a very thin sheet metal made of stainless steel with a thickness of less than 200 $\mathsf{\mu }$m. Given these situations, a set of simplifying assumptions is outlined in the following similar to the reference model introduced by [10,13].

• The material is isotropic and homogeneous with constant thermo-physical characteristics.
• The input laser follows a Gaussian beam distribution and is considered to be a surface heat flux.
• The laser beam is absorbed on the surface and emitted vertically downward. The absorptivity coefficient is taken to be constant and describes the absorbed power over the surface projected in the direction of the laser beam.
• The laser beam radius is constant (=beam waist).
• The melt flow, plasma formation, and recoil pressure are not considered in the model.
• The molten material point is removed as soon as it reaches the melting temperature (without considering the latent heat) and does not show hydrodynamic behavior.
• The removed material is transparent and does not absorb the laser energy.

The corollary to these simplifications is to allow for a single-phase thermal modeling of the laser drilling process.

2.2. Governing Equations

The thermal field is governed by the Fourier heat equation for isotropic solids expressed in a 3D Cartesian coordinate system as:

$$\rho c_p\frac{\partial T}{\partial t}=k{\nabla }^{2}T=k\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}\right)$$

(1)

where $\rho$ is the density, $c_p$ specific heat capacity, T temperature, and k heat conductivity. The boundary conditions associated with Equation (1) are the surface heat flux terms defined as:

• Heat source $q_s$ due to the absorbed laser beam.
• Heat loss $q_l$ due to convection and radiation.

These terms are given by:

$$\begin{array}{cc}\hfill q_s=& {\displaystyle n_z\frac{a_L{P}_L}{\pi r_b^2}exp\left[-2{\left(\frac{r}{r_b}\right)}^2\right]}\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill q_l=& {\displaystyle \underset{\mathrm{convection}}{\underbrace{h_c(T_s-T_{\infty })}}+\underset{\mathrm{radiation}}{\underbrace{ϵ\sigma (T_s^4-T_{\infty }^4)}}}\hfill \end{array}$$

(3)

in which $n_z$ is the direction cosine of the normal vector outward from the laser-irradiated surface, $a_L$ the laser absorption coefficient ($0<a_L<1$), $P_L$ the laser peak power, r the radial distance from the beam center, and $r_b$ the laser beam radius (i.e., in the waist). In Equation (3), $h_c$ is the heat convection coefficient, $T_s$ and $T_{\infty }$ the surface and background temperatures, respectively, $ϵ$ the emissivity factor, and $\sigma$ the Stefan–Boltzmann constant.

3. Numerical Modeling

An overview of the numerical methods to be used for the thermal model described in Section 2 is the central focus of this section. For more clarity, the problem is decomposed into the following subsets.

• Discretization of the Laplace operator in Equation (1).
• Treatment of boundary conditions in Equations (2) and (3).
• Time integration of Equation (1).

Details regarding the temporal discretization of an ODE in meshfree models are not provided in this paper. It suffices to point out that an explicit leapfrog stepping according to [29] is used for the time integration unless otherwise stated. It is widely accepted in explicit algorithms to impose a maximum step size $\Delta t$ for the sake of solution stability. In this work, the step size is calculated from the Courant–Friedrichs–Lewy (CFL) criterion by restricting the maximum permissible value of $\Delta t$ to obey:

$$\Delta t\le \chi \frac{{(\Delta x)}^2}{\alpha }$$

(4)

where $\Delta x$ is the characteristic distance between particles, $\alpha =k/\rho c_p$ the thermal diffusivity, and $\chi$ a scale factor considered to be 0.25 in this paper. Discussions regarding the first and second parts of the problem decomposition are presented with more detail in the next two subsections.

3.1. Meshfree Schemes for Laplacian Approximation

The present framework requires no computational mesh for spatial discretizations and is built on the notion of SPH techniques. In principle, SPH is an interpolation method based on the position of Lagrangian particles and the use of a smooth weighting function called the kernel. Following the original derivation by [30,31], the numerical procedure of SPH for an arbitrary field $f\left(\underset{¯}{r}\right)$ at particle i consists of two approximations:

$$\begin{array}{cc}\hfill f({\underset{¯}{r}}_i)\approx \langle f\left({\underset{¯}{r}}_i\right)\rangle ={\int }_{\mathrm{\Omega }}f\left(\underset{¯}{r}\right)W\left(\right|{\underset{¯}{r}}_i-\underset{¯}{r}|,h)dV& \to \mathrm{kernel}\mathrm{approximation}\left(\mathrm{interpolation}\right)\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \langle f\left({\underset{¯}{r}}_i\right)\rangle \approx \sum _{j}f\left({\underset{¯}{r}}_j\right)W\left(\right|{\underset{¯}{r}}_i-{\underset{¯}{r}}_j|,h)V_j& \to \mathrm{particle}\mathrm{approximation}\left(\mathrm{discretization}\right)\hfill \end{array}$$

(6)

where $\mathrm{\Omega }$ is the computational domain, $dV$ a differential volume element, and W a kernel function that depends only on two variables: the particle distance ${\underset{¯}{r}}_{ij}=|{\underset{¯}{r}}_{ij}|=|{\underset{¯}{r}}_i-{\underset{¯}{r}}_j|$ and the smoothing length h. A variety of frequently-used kernels in SPH models were listed in [29,32], from which the following functions were chosen for the numerical examples of this paper.

• Cubic B-spline:

  $$W_h\left(\underset{¯}{r}\right)=n_c\times \left\{\begin{array}{cc}{\displaystyle 1-\frac{3}{2}{q}^2+\frac{3}{4}{q}^3,}\hfill & 0\le q<1\hfill \\ {\displaystyle \frac{1}{4}{(2-q)}^3,}\hfill & 1\le q\le 2\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

  (7)
• Wendland quintic:

  $$W_h\left(\underset{¯}{r}\right)=n_c\times \left\{\begin{array}{cc}{\displaystyle {(1-\frac{q}{2})}^4(2q+1),}\hfill & q<2\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

  (8)

where r is the position vector, h the smoothing length, which controls the amount of mollification (smoothing), and $q=\left|\underset{¯}{r}\right|/h$. The normalization constant for 3D space is $n_c=1/\pi h^3$ and $n_c=21/16\pi h^3$ for the cubic B-spline and Wendland quintic kernels, respectively.

A major novelty of this article lies in the first application of these methods to the (3D) laser drilling problem, while the 3D implementation of ICSPM is presented for the first time. Without delving deeply into their derivations, the six variants listed in Table 1 of the introductory section are reviewed in consistent notation and a self-contained fashion.

3.1.1. Brookshaw’s Approximation (BSH)

The method presented by Brookshaw [22] is considered as one of the most popular options in SPH implementations. It benefits from a simple correction to the original SPH form and is expressed by:

$$\begin{array}{c}\hfill {\langle {\nabla }^2{T}_i\rangle }^{\mathrm{BSH}}\approx \sum _{j}2\left[\left(\frac{T_i-T_j}{|r_{ij}|}\right){\underset{¯}{e}}_{ij}·\nabla W_{ij}\right]V_j\end{array}$$

(9)

where ${\underset{¯}{e}}_{ij}={\underset{¯}{r}}_{ij}/\left|r_{ij}\right|$ is a unit vector in the inter-particle direction. BSH guarantees the conservation of energy to the precision of the time-stepping algorithm, which is an important feature of this scheme.

3.1.2. Particle Strength Exchange (PSE)

Introduced by [23,33] for the Laplacian approximation of advection-diffusion problems, PSE is regarded as a strong numerical tool for solving diffusion problems. Instead of computing the derivatives of the same SPH kernel, PSE constructs different kernels for different derivatives independently. A good review of this method can be found in [34,35]. Discretization of the three-dimensional Laplacian of T using a second-order PSE kernel reads:

$${\langle {\nabla }^2{T}_i\rangle }^{\mathrm{PSE}}\approx \sum _j(T_j-T_i)\frac{1}{h^3}\left[\frac{4}{\pi \sqrt{\pi }{h}^2}exp\left(-\left(\right|{\underset{¯}{r}}_{ij}{|/h)}^2\right)\right]V_j$$

(10)

where the pre-factor $1/h^3$ appears due to the scaling property of PSE kernels in 3D geometries. A more exhaustive list of other PSE kernels for various differential operators and different dimensions can be found in the Appendix of [34]. A comparison between Equations (9) and (10) reveals that the PSE approximation also offers the conservation of energy up to the precision of the time-stepping algorithm. It is understood from Equations (9) and (10) that the Laplacian results using BSH and PSE are not generally consist because of their kernel deficiency near/on the boundaries. The corrective schemes presented in the next subsections are employed to address this issue.

3.1.3. Reproducing Kernel Particle Method (RKPM)

RKPM is a meshfree technique, which, in theory, can be derived for a favorable order of accuracy and consistency by incorporating different correction functions. In 1995, Liu et al. [24] originally applied this method in their seminal work to advection-diffusion problems. The constructive idea of RKPM is to enhance SPH kernels by applying a correction function $\psi$ for recovering the consistency. Iterating the same procedure described in [35] for 3D implementations leads to a first-order consistent RKPM approximation for second derivatives, and subsequently, for the Laplace operator. If the approximation of ${\partial }^2/\partial x^2$ for T is desired, the mathematical expression arrives at:

$$\begin{array}{c}\hfill {\displaystyle {〈\frac{{\partial }^2{T}_i}{\partial x^2}〉}^{\mathrm{RKPM}}\approx \sum _j{T}_j\left[\left(\frac{{\partial }^2{\psi }_{ij}}{\partial x^2}\right)W_{ij}+2\left(\frac{\partial {\psi }_{ij}}{\partial x}\right)·\left(\frac{\partial W_{ij}}{\partial x}\right)+{\psi }_{ij}\left(\frac{{\partial }^2{W}_{ij}}{\partial x^2}\right)\right]V_j}\end{array}$$

(11)

where $\psi$ is a linear correction function whose unknown coefficients are determined by enforcing a set of moment conditions and solving a system of linear equations. This procedure is repeated for the other spatial coordinates to find the Laplacian value by summing three individual second derivatives (in 3D). As pointed out in [35], it is quite popular to utilize a symmetric or anti-symmetric form of Equation (11) by replacing $T_j$ with $[T_j\pm T_i]$.

3.1.4. Corrective Smoothed Particle Method (CSPM)

CSPM was initially derived in [25], but presented later by Chen and Beraun [36] in a more generalized framework. The constructive idea of CSPM is to reduce errors in the first and second derivatives computed with SPH by correcting the gradient operators and manipulating the Taylor series expansion. For the calculation of the second derivatives in CSPM, the values of the first derivatives must be known beforehand. The CSPM discretization of the second derivatives involves multiple matrix operations and a system of linear equations for each particle at each time step. It can be written as:

$$\begin{array}{cc}\hfill {〈\nabla \otimes \nabla T_i〉}^{\mathrm{CSPM}}& \approx {\underset{\approx }{\mathit{B}}}_i^{-1}\left[\sum _j(T_j-T_i)\nabla \otimes \nabla W_{ij}{V}_j-\sum _j\left({\underset{¯}{r}}_{ji}·〈\nabla T_i〉\right)\nabla \otimes \nabla W_{ij}{V}_j\right]\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill 〈\nabla T_i〉& \approx {\underset{~}{\mathit{A}}}_i^{-1}\left[\sum _j(T_j-T_i)\nabla W_{ij}{V}_j\right]\hfill \end{array}$$

(13)

whereby the fourth-order tensor $\underset{\approx }{\mathit{B}}$ is defined for particle i:

$$\begin{array}{cc}\hfill {\underset{\approx }{\mathit{B}}}_i& =\sum _j\left({\underset{~}{\mathit{C}}}_{ij}\otimes \nabla \otimes \nabla W_{ij}\right)V_j\hfill \end{array}$$

(14)

in which the expanded definition of matrix $\underset{~}{\mathit{C}}$ is found in [35,36], and the renormalization matrix $\underset{~}{\mathit{A}}$ for particle i is computed by:

$$\begin{array}{c}\hfill {\underset{~}{\mathit{A}}}_i=\sum _j\left({\underset{¯}{r}}_{ji}\otimes \nabla W_{ij}\right)V_j\end{array}$$

(15)

It should be noted that $\langle \nabla \otimes \nabla T\rangle$ in Equation (12) indicates the Hessian matrix of T, meaning that the numerical value of ${\nabla }^{2}T$ still needs to be calculated by the sum of elements on the main diagonal. Chen et al. [25,37] showed that CSPM is a first-order consistent method that enhances the solution accuracy, not only near the boundary, but also inside the domain.

3.1.5. Improved CSPM

More recently, Korzilius et al. [27] showed that a second-order consistent discretization of the Laplace operator can be obtained by incorporating a discrete first derivative, which is second-order consistent. They included one more term in the corresponding Taylor series of $\nabla T$ in Equation (13) and derived a new fourth-order renormalization tensor $\underset{\approx }{\mathit{D}}$ for the second derivatives:

$${〈\nabla \otimes \nabla T_i〉}^{\mathrm{ICSPM}}\approx {\underset{\approx }{\mathit{D}}}_i^{-1}\left[\sum _j(T_j-T_i)\nabla \otimes \nabla W_{ij}{V}_j-\sum _j\left({\underset{¯}{r}}_{ji}·\langle \nabla T_i\rangle \right)\nabla \otimes \nabla W_{ij}{V}_j\right]$$

(16)

where $\langle \nabla T_i\rangle$ was defined in Equation (13) and $\underset{\approx }{\mathit{D}}$ is calculated from:

$${\underset{\approx }{\mathit{D}}}_i={\underset{\approx }{\mathit{B}}}_i-\left[\left(\sum _j({\underset{¯}{r}}_{ji}\otimes \nabla \otimes \nabla W_{ij})V_j\right)·{\underset{~}{\mathit{A}}}_i·\left(\sum _j({\underset{~}{\mathit{C}}}_{ij}\otimes \nabla W_{ij})V_j\right)\right]$$

(17)

Just like CSPM, the Laplacian value using ICSPM needs to be calculated by the sum of the elements on the main diagonal of Equation (16). The ICSPM approximation gives a second-order consistent solution [27], while its extra computational effort compared to CSPM is negligible. An adaptation of ICSPM for thermal modeling of 2D multi-phase flows was presented by [38].

3.1.6. Fatehi’s Meshfree Scheme (FMFS)

Originally derived by Fatehi et al. [26] and recently applied to 3D heat conduction problems in [35], FMFS is another modern SPH scheme that provides a second-order consistent approximation of the Laplace operator. It computes $\langle {\nabla }^2\rangle$ of function T via:

$${〈{\nabla }^2{T}_i〉}^{\mathrm{FMFS}}\approx {\underset{~}{\mathit{G}}}_i:\sum _{j}2\left[\left(\frac{(T_j-T_i)}{|{\underset{¯}{r}}_{ij}|}-{\underset{¯}{e}}_{ij}·\sum _j{V}_j(T_j-T_i){\underset{~}{\mathit{A}}}_i^{-1}·\nabla W_{ij}\right){\underset{¯}{e}}_{ij}\otimes \nabla W_{ij}\right]V_j$$

(18)

by introducing a novel second-order renormalization tensor $\underset{~}{\mathit{G}}$ computed from:

$${\underset{~}{\mathit{G}}}_i:\left[{\displaystyle {\underset{~}{\underset{¯}{\mathit{E}}}}_i·{\underset{~}{\mathit{A}}}_i^{-1}·{\underset{~}{\underset{¯}{\mathit{F}}}}_i+\sum _j\left({\underset{¯}{r}}_{ij}\otimes {\underset{¯}{e}}_{ij}\otimes {\underset{¯}{e}}_{ij}\otimes \nabla W_{ij}\right)V_j}\right]=-\underset{~}{\mathit{I}}$$

(19)

in which “:” denotes the double contraction of two second-order tensors. The CSPM renormalization tensor $\underset{~}{\mathit{A}}$ was defined in Equation (15), and ${\underset{¯}{e}}_{ij}=({\underset{¯}{x}}_{ij}+{\underset{¯}{y}}_{ij}+{\underset{¯}{z}}_{ij})/r_{ij}$ is the normalized or unit vector. The terms $\underset{~}{\mathit{E}}$ and $\underset{~}{\mathit{F}}$ in Equation (19) are third-order tensors given by:   

$$\begin{array}{cc}\hfill {\underset{~}{\underset{¯}{\mathit{E}}}}_i=& \sum _j\left({\underset{¯}{e}}_{ij}\otimes {\underset{¯}{e}}_{ij}\otimes \nabla W_{ij}\right)V_j\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\underset{~}{\underset{¯}{\mathit{E}}}}_i=& \sum _j\left({\underset{¯}{r}}_{ij}\otimes {\underset{¯}{r}}_{ij}\otimes \nabla W_{ij}\right)V_j\hfill \end{array}$$

(21)

which are sensitive only to the location of particles. In order to minimize the computational effort of this expensive method, it is beneficial to maximize the amount of pre-computation intensity. Because the particles’ location is the only ingredient that constitutes the correction terms and renormalization tensors, the pre-computation block (which appears in corrected SPH methods) needs to be called only if any changes occur to the particle arrangement. This statement applies to other corrected SPH methods as well.

3.2. Proposed Approach for Boundary Issues

A simple and efficient strategy is proposed for the identification of surface particles, which can be used for thermal modeling applications without material deformations. The rationale behind this approach is twofold: (1) to calculate the unit normal vector at the laser-irradiated surface; (2) to impose the boundary conditions defined in Equations (2) and (3) on surface particles at each time-step.

Since the source term $q_l$ is a surface heat flux, the direction cosine of the normal vector outward from the laser-irradiated surface, $n_z$ in Equation (2), needs to be known. Calculation of the unit normal vector in this work follows the color function idea used for the surface tension term in CFD applications [16,38]. For particle i:

$${\widehat{\underset{¯}{n}}}_i=\frac{{\underset{¯}{n}}_i}{|{\underset{¯}{n}}_i|}\approx \frac{\nabla c_i}{|\nabla c_i|}$$

(22)

wherein $\widehat{\underset{¯}{n}}$ is the unit normal vector and $\nabla c$ the gradient of the color function calculated via a symmetric SPH approximation:

$$\nabla c_i\approx \frac{1}{V_i}\sum _j\left(V_i^2+V_j^2\right)\frac{c_i^i+c_i^j}{2}\nabla W_{ij}$$

(23)

with V being the integration volume of particles. This is in fact a density-weighted color function [16], in which the index number $c_i^j$ between two different phases (or colors) is computed from:

$$c_i^j=\left\{\begin{array}{cc}{\displaystyle \frac{2{\rho }_i}{{\rho }_i+{\rho }_j}}\hfill & \mathrm{if}\mathrm{particle}i\mathrm{and}j\mathrm{belong}\mathrm{to}\mathrm{different}\mathrm{phases}/\mathrm{labels}\hfill \\ 0\hfill & \mathrm{if}\mathrm{particle}i\mathrm{and}j\mathrm{belong}\mathrm{to}\mathrm{the}\mathrm{same}\mathrm{phase}/\mathrm{label}\hfill \end{array}\right.$$

(24)

Moreover, one notices that $q_l$ should be applied only to the laser-irradiated particles. A straightforward approach is suggested for detecting these surface particles during the simulation. It performs locally and makes use of a state flag for each particle. In summary, the order of operations is given in the following four steps.

• Particle i is marked molten as soon as its temperature reaches the melting point.
• The neighbor particles adjacent to particle i are marked with a candidate flag.
• All candidate particles that are neither molten nor already surface particles become new surface particles.
• All molten particles are deleted.

Clearly, this strategy is valid (and efficient) as long as particles are fixed in space and the array of particle positions remains unchanged in memory, ensuring a fast identification of the neighboring particles. Figure 2 illustrates how the proposed approach works at a given time-step. Due to symmetry and for the sake of better visualization, only half of the workpiece is displayed. It is readily apparent that this laser drilling problem has a high symmetry around the irradiation axis, hence is a 2D axisymmetric model. Nevertheless, the present simulation was run in 3D to verify the computer implementation and efficiency.

4. Results and Discussion

Before simulating the laser drilling model in this section, a validation benchmark on a 3D transient heat conduction problem was conducted to ensure the correct implementation of the presented methods. As for the computer software, the open-source code “Thermal-IWF” written in C99 and published by [35] is utilized. To download this solver, please refer to the https://github.com/mroethli/thermal_iwf repository. All computations were performed on an Intel${}^{\circledR }$ Core™ i5-4690 with four processors at 3.50 GHz. Most of the images were created with ParaView [39]—a powerful visualization tool for large datasets.

4.1. Validation Benchmark

The test case was related to the time-dependence of temperature fields in a 3D object and examined the accuracy and performance of the Laplacian models presented in Section 3.1. If the initial condition of Equation (1) follows a thermal pulse at the origin of the domain $(x_o,y_o,z_o)$ without Dirichlet boundary conditions, the time-dependent analytical solution of this problem at time t can be computed from:

$${\displaystyle T=\frac{1}{\left(4\pi \alpha t\right)\sqrt{\left(4\pi \alpha t\right)}}exp\left(\frac{-{\left|r\right|}^2}{4\alpha t}\right)}$$

(25)

in which $\alpha$ is the thermal diffusivity assumed to be one in this model, $\left|r\right|$ the radial distance from the origin, and T the temperature value at time t generated by an instantaneous source of heat at the origin point $P_o(x_o,y_o,z_o)$ at time $t=t_0$. The time span for the transition of heat in this mode was considered to start from $t_0=0.1$ and end at $t_f=0.3$. The infinite boundaries were herein limited to a cube of $x\in \left[-5,+5\right],y\in \left[-5,+5\right],z\in \left[-5,+5\right]$, as shown by the black outlined boxes in Figure 3. This domain size was deemed sufficiently large because the boundary temperature due to the heat impulse in the origin at $t_f$ would be about ${10}^{-19}$, orders of magnitude less than the error due to the methods and, in fact, orders of magnitude smaller than the machine epsilon.

The particle spacing (i.e., discretization size) was chosen from $\Delta x=\{1,0.5,0.25\}$ to evaluate the errors and draw the convergence plots. All simulations were carried out on both regular and random particle arrangements to compare the capabilities of the methods in handling disordered situations. The randomness of particles was limited to the maximum of 50% inhomogeneity for all particles except from the edge points. The edge particles were always distributed with a regular spacing (see Figure 3).

The numerical results were compared with the exact solution at $t_f=0.3$, and the ${\mathcal{L}}_2$ norm of the error is reported. The convergence plots shown in Figure 4 permit several interesting observations. Firstly, it can be seen that almost all methods converged with a second-order slope in the regular spacing (left graph in Figure 4). Secondly, CSPM and ICSPM delivered identical results for the uniform case. This was quite unsurprising because the main advantage of ICSPM over CSPM is its ability to provide higher-order consistency by improving the boundary deficiency, but the boundary effect in this example was negligible. When particles were not uniformly spaced, however, a marginal difference between the CSPM and ICSPM results was noticed. Thirdly, FMFS demonstrated the lowest error and outperformed the other Laplacian models in almost all cases. This was due to the higher order terms included in its truncation error. A particularly interesting conclusion from this investigation was the nearly identical behavior of SPH Brookshaw (the red graphs in Figure 4) and FMFS in discretizing the Laplacian term for uniform cases. The reason lied in the derivation process of FMFS, which originated from a Brookshaw approximation. Please consider [35] for more details in this regard.

4.2. Laser Drilling Simulation

The case study pertained to a sheet metal made of stainless steel (SS 316L) with a width of $l_x=100$ $\mathsf{\mu }$m, length of $l_y=200$$\mathsf{\mu }$m, and thickness of $l_z=150$$\mathsf{\mu }$m. Consistent with the initial configuration of the 3D model used by Muhammad et al. [10], the same material and laser properties were taken into account. Listed in

Table 2. Parameters and material properties used according to [10,13].

Property	Unit	Value
Density	kg m${}^{-3}$	7950
Specific heat capacity	J kg${}^{-1}$ K${}^{-1}$	470
Heat conductivity	W m${}^{-1}$ K${}^{-1}$	20
Heat convection coefficient	W m${}^{-2}$ K${}^{-1}$	20
Stefan–Boltzmann constant	W m${}^{-2}$ K${}^{-4}$	5.67×10${}^{-8}$
Emissivity factor	–	0.4
Initial temperature	K	300
Melting temperature	K	1723
Laser power	W	100
Laser absorptivity	–	0.4
Laser beam radius	$\mathsf{\mu }$m	25

is a summary of the parameters used for this simulation. The absorptivity coefficient $a_L$ was chosen to be $40\%$ according to a very relevant investigation performed by [13] on the sensitivity of MLPG results to this parameter.

Muhammad et al. [10] were the first who considered this experimental setup for SPH modeling and showed that a pulse of 150 $\mathsf{\mu }$s could fully penetrate the metallic sheet with a thickness of 150 $\mathsf{\mu }$m. Three images adapted from the laser cutting experiments of [10] are presented in Figure 5 for three different pulse durations. As a result of this observation, the present thermal simulation was terminated at $t=150$ $\mathsf{\mu }$s to investigate the computed temperature distribution, removed volume, and penetration depth.

Firstly, a parametric study of this example was carried out in 2 spatial resolutions using 2 kernel functions and 3 different smoothing lengths, for all 6 methods presented in Section 3.1. This combination led to a massive load of simulations including 72 configuration setups. The numerical results are inserted into

Table 3. Batch 1: Coarse discretization using a total of 27k particles with the cubic B-spline kernel.

Method	Data	$\mathit{h}=\Delta \mathit{x}$	$\mathit{h}=1.25\Delta \mathit{x}$	$\mathit{h}=1.5\Delta \mathit{x}$
BSH	Rem.vol.(%)	63.41	74.36	79.37
	Full pen.	✓	✓	✓
	Runtime (s)	32	60	102
PSE	Rem. vol (%)	72.83	78.13	83.4
	Full pen.	✓	✓	✓
	Runtime (s)	77	137	237
RKPM	Rem. vol. (%)	63.28	63.53	78.18
	Full pen.	✓	✓	✓
	Runtime (s)	28	53	86
CSPM	Rem. vol. (%)	74.99	83.82	–
	Full pen.	✓	✓	–
	Runtime (s)	73	130	200
ICSPM	Rem. vol. (%)	74.99	83.02	–
	Full pen.	✓	✓	–
	Runtime (s)	73	131	200
FMFS	Rem. vol. (%)	80.77	80.68	84.76
	Full pen.	✓	✓	✓
	Runtime (s)	38	74	123

,

Table 4. Batch 2: Coarse discretization using a total of 27k particles with the Wendland quintic kernel.

Method	Data	$\mathit{h}=\Delta \mathit{x}$	$\mathit{h}=1.25\Delta \mathit{x}$	$\mathit{h}=1.5\Delta \mathit{x}$
BSH	Rem. vol. (%)	65.95	72.66	78.56
	Full pen.	✓	✓	✓
	Runtime (s)	32	61	101
PSE	Rem. vol (%)	72.83	78.13	83.4
	Full pen.	✓	✓	✓
	Runtime (s)	77	137	237
RKPM	Rem. vol. (%)	78.47	76.44	80.47
	Full pen.	✓	✓	✓
	Runtime (s)	29	53	88
CSPM	Rem. vol. (%)	76.52	77.63	87.85
	Full pen.	✓	✓	✓
	Runtime (s)	73	129	205
ICSPM	Rem. vol. (%)	76.52	77.67	89.51
	Full pen.	✓	✓	✓
	Runtime (s)	73	129	205
FMFS	Rem. vol. (%)	80.47	81.44	85.39
	Full pen.	✓	✓	✓
	Runtime (s)	38	74	122

,

Table 5. Batch 3: Fine discretization using a total of 203k particles with the cubic B-spline kernel.

Method	Data	$\mathit{h}=\Delta \mathit{x}$	$\mathit{h}=1.25\Delta \mathit{x}$	$\mathit{h}=1.5\Delta \mathit{x}$
BSH	Rem. vol. (%)	40.17	39.35	68.66
	Full pen.	–	–	✓
	Runtime (s)	671	1234	2095
PSE	Rem. vol (%)	60.21	74.92	86.73
	Full pen.	✓	✓	✓
	Runtime (s)	2405	4489	7628
RKPM	Rem. vol. (%)	32.61	33.42	–
	Full pen.	–	–	–
	Runtime (s)	1134	2167	3536
CSPM	Rem. vol. (%)	71.50	–	–
	Full pen.	✓	–	–
	Runtime (s)	2589	4541	5357
ICSPM	Rem. vol. (%)	71.50	–	–
	Full pen.	✓	–	–
	Runtime (s)	2592	4515	5609
FMFS	Rem. vol. (%)	87.74	89.63	90.69
	Full pen.	✓	✓	✓
	Runtime (s)	1358	2653	4433

and

Table 6. Batch 4: Fine discretization using a total of 203k particles with the Wendland quintic kernel.

Method	Data	$\mathit{h}=\Delta \mathit{x}$	$\mathit{h}=1.25\Delta \mathit{x}$	$\mathit{h}=1.5\Delta \mathit{x}$
BSH	Rem. vol. (%)	41.48	39.08	63.13
	Full pen.	–	–
	Runtime (s)	660	1268	2115
PSE	Rem. vol (%)	60.21	74.92	86.73
	Full pen.	✓	✓	✓
	Runtime (s)	2405	4489	7628
RKPM	Rem. vol. (%)	87.03	46.69	48.91
	Full pen.	✓	✓	✓
	Runtime (s)	1106	2015	3364
CSPM	Rem. vol. (%)	76.51	69.01	97.55
	Full pen.	✓	✓	✓
	Runtime (s)	2022	4666	7630
ICSPM	Rem. vol. (%)	76.51	63.43	–
	Full pen.	✓	✓	–
	Runtime (s)	2028	4673	7591
FMFS	Rem. vol. (%)	87.74	89.42	90.97
	Full pen.	✓	✓	✓
	Runtime (s)	1213	2819	4534

, where quantities of three different kinds are reported for assessment:

• The percentage of the volume removed from the workpiece after $\Delta t=150\mathsf{\mu }$s if the total volume to be removed was conceptually equal to a cylinder of radius $r_b=25\mathsf{\mu }$m and height $l_z=150$$\mathsf{\mu }$m.
• Whether or not a full penetration was observed.
• The measured CPU runtime of the simulation.

The parameter $\Delta x$ in these four tables denotes the characteristic distance between particles in a uniform domain. Furthermore, the appearance of “—” corresponds to a scenario where the numerical result was unreliable or invalid. The source of these errors may stem from different issues such as the surface detection algorithm, boundary conditions, insufficient number of neighbor particles, and inconsistent kernels. An elaboration of how to resolve such issues is irrelevant to the scope of this paper and thus not provided. For further discussions, the readers are referred to textbooks and review articles such as [40,41]. Note that $\langle {\nabla }^{2}T\rangle$ at surface particles was always computed by PSE even if other corrective schemes (RKPM, CSPM, ICSPM, or FMFS) were employed in the code. This new workaround ensured the stability of the solution without requiring exterior dummy particles for an invertible correction matrix where boundaries were not well defined geometrically (e.g., the laser-irradiated surface inside the popped hole).

As can be seen in these tables, CSPM and ICSPM showed the least stable behavior among the other schemes. This issue was due to the excessive number of inverse matrix calculations and also interrelated with the proposed boundary treatment approach. Unsurprisingly, the choice of the kernel function had a major impact on the stability of the simulation results. A comparison between Table 3 and Table 4 in the case of CSPM/ICSPM supported the argument that the cubic B-spline kernel performed weaker than its contestant at higher h. As a result, the Wendland quintic kernel function was selected for final simulations at a high resolution.

As far as the computational cost was concerned, attention must be paid when comparing the runtime of different schemes in purely thermal simulations. In all corrective methods (RKPM, CSPM, ICSPM, and FMFS), a significant proportion of the computational cost was attributed to the computation of their correction matrices. However, the evaluation of these matrices here was not updated within the main time loop, and a remarkable amount of computational time could be saved because the particle position remained unchanged during a purely thermal simulation. Interestingly though, the PSE scheme was indifferent to the movement of particles. It considered an exponential kernel according to Equation (10), no matter if the particle position had been changed or not. Therefore, any conclusive judgment about the computational performance of these six methods was reckoned to be inaccurate and was probably unfair.

After having completed the parametric study, the optimum setup for each method could be configured. In Figure 6, the removed volume predicted by this 3D solver is compared to the corresponding value taken from a 2D meshless model developed by Abidou et al. [13]. The bar-chart shows that the highest volumes were computed by FMFS on account of its higher order terms involving a larger number of particles for each interpolation.

Due to the sufficient simplifications considered in this test model, any validation or comparison versus experimental measurements would be questionable at best. Nevertheless, it was still worthwhile to include a demonstration of the present result against other references. Simulated penetration depths computed by this code are plotted in Figure 7 for all six methods at two different resolutions. The numerical values and overall trend were found acceptable compared to the experimental data provided by [10] and the numerical solutions of [13]. It is important to mention that the numerical simulation presented by [13] was also performed by taking a set of assumptions similar to what was outlined in this paper.

The experimental data points in Figure 8 were the penetration depth recorded for each pulse duration. In this figure, the presented methods were alternatively used to predict the penetration in time. A more or less linear behavior was noticed in all numerical graphs. According to the experimental evidence, however, the laser absorptivity coefficient at lower temperatures was less than its values at higher temperatures [10,12,13]. In other words, at the beginning of the laser drilling process, this coefficient had a lower value compared to the later stages, when it abruptly increased due to the rapid increase in both temperature and depth. Another reason for this behavior was that the laser beam was also defocused in lower parts of the hole, which indicated lower intensity and, therefore, lower heat input and surface temperature. In fact, the penetration speed was slow in the beginning as no multiple reflections occurred, and the laser due to the orthogonal irradiation had a low absorptivity. The speed increased due to multiple reflections and much better absorptivity on the steep walls. These were all the effects that were ignored by the simplifying assumptions in this paper and must be taken into consideration for future development.

Next, we compared the computational performance of the presented methods by reporting their runtime at three different resolutions of 27 k, 203 k, and $1.6$ mtotal particles. Both linear and log–log graphs are plotted in Figure 9, demonstrating the highest cost of computation for CSPM/ICSPM and the lowest for BSH. Furthermore, an exemplary snapshot of each resolution is represented in Figure 10 to illustrate the impact of (spatial) discretization size on the quality of the results.

Lastly, images of the temperature distribution are provided in high resolution to gain more insights into the process. Figure 11 contains six frames of the simulation taken at six consecutive time instants, where PSE was employed for the discretization of the Laplacian operator. As a verification that our results were correct, the pictures in Figure 11 could be compared with those generated by [10,13] for a very similar case. These snapshots also served to represent how the laser-induced heat was transferred into the metallic workpiece, leading to the formation of the drilled hole. A summary of this scenario, which was captured by the present implementation rather nicely, follows. Initially, top surface particles absorbed the thermal energy provided by the laser beam. Due to the heat conduction, material particles were heated up until they reached the melting point and consequently become inactive in the simulation. These molten particles (invisible in Figure 11 because they were inactive in the model) were removed from the popped hole, and the laser beam was applied onto a newly generated surface (i.e., the orange particles in Figure 2). Therefore, the evolutionary boundary conditions during this procedure were effectively handled; the drilled hole was progressively developed; and the drilling penetration was obtained.

5. Conclusions

This paper developed a 3D numerical simulation framework that could return quantitative and qualitative insights into the thermal modeling of laser drilling processes. Six meshfree schemes with different orders of consistency were first cross-compared in a convergence study and then applied for the first time to the 3D thermal modeling of a laser drilling problem. While these methods shared some similarities in the quadratic rate of convergence for the approximation of the Laplace operator (see Figure 4), their computational cost and stability were fundamentally different in the laser drilling simulation (see Table 3, Table 4, Table 5 and Table 6 and Figure 9, for example). It was shown in Figure 6 that the total removed volume calculations compared well with those reported by [13] using the MPLG collocation method. Furthermore, the comparison of available experimental data from [10] and the penetration depth computed by this solver were found to be within an acceptable range considering the simplified modeling procedure. Nonetheless, model refinements and further experimental studies are envisaged to be necessary for better evaluation.

Overall, the computer code performed efficiently, and the implemented methods could provide useful insights into the thermal field and fundamental nature of a laser drilling process. However, an immediate enhancement to this work would adopt a more sophisticated approach for handling general boundary conditions where particles can freely move. Moreover, it is of utmost importance to include more physical phenomena in the description to increase the degree of realism required for a more realistic model. For instance, accounting for the change of laser absorptivity, modeling of the phase change, latent heat, the Marangoni effect, assist gas, molten particles, and recoil pressure could be added to future developments.