# An Incompressible Smoothed Particle Hydrodynamics (ISPH) Model of Direct Laser Interference Patterning

^{1}

^{2}

^{*}

Next Article in Journal

Previous Article in Journal

Institute of Manufacturing Technology, Technische Universität Dresden, P.O. Box, 01062 Dresden, Germany

Fraunhofer Institute for Material and Beam Technology IWS, Winterbergstraße 28, 01277 Dresden, Germany

Author to whom correspondence should be addressed.

Received: 19 December 2019 / Revised: 20 January 2020 / Accepted: 27 January 2020 / Published: 30 January 2020

(This article belongs to the Section Computational Engineering)

Functional surfaces characterised by periodic microstructures are sought in numerous technological applications. Direct laser interference patterning (DLIP) is a technique that allows the fabrication of microscopic periodic features on different materials, e.g., metals. The mechanisms effective during nanosecond pulsed DLIP of metal surfaces are not yet fully understood. In the present investigation, the heat transfer and fluid flow occurring in the metal substrate during the DLIP process are simulated using a smoothed particle hydrodynamics (SPH) methodology. The melt pool convection, driven by surface tension gradients constituting shear stresses according to the Marangoni boundary condition, is solved by an incompressible SPH (ISPH) method. The DLIP simulations reveal a distinct behaviour of the considered substrate materials stainless steel and high-purity aluminium. In particular, the aluminium substrate exhibits a considerably deeper melt pool and remarkable velocity magnitudes of the thermocapillary flow during the patterning process. On the other hand, convection is less pronounced in the processing of stainless steel, whereas the surface temperature is consistently higher. Marangoni convection is therefore a conceivable effective mechanism in the structuring of aluminium at moderate fluences. The different character of the melt pool flow during DLIP of stainless steel and aluminium is confirmed by experimental observations.

Microscopic periodic features provide surfaces with superior functionalities. In biological and medical applications, repetitive surface textures improve the biocompatibility of bone implants [1], guide directional cell growth [2] and inhibit bacterial adhesion and biofilm formation [3]. Further advantages of periodic structured surfaces include enhanced light absorption [4], reduced friction [5,6] and anisotropic wetting [7]. These topographies are increasingly considered for the manufacture of functional surfaces, e.g., in biomedical and marine engineering, tribology, optics and aeronautics.

Direct laser interference patterning (DLIP) is a novel method that allows the production of periodic surface structures with feature sizes in the submicron to micron range in a single processing step. For this purpose, the primary beam of a pulsed laser with wavelength $\lambda $ is split into two or more partial beams. The periodic intensity distribution due to the interfering coherent beams is employed to treat a surface situated in the interference volume. Here, the interference of two beams is considered, resulting in a sinusoidal energy density distribution [8]
where ${\Phi}_{0}$ is the fluence of each beam and $\theta $ is the intersecting angle between the beams. Accordingly, two-beam interference patterning generates line-like surface structures with a spatial periodicity
the minimum achievable periodic distance being half of the used wavelength, i.e., $\Lambda =\lambda /2$ for $\theta =\mathsf{\pi}$. The fabrication of repetitive microstructures by means of DLIP, using nanosecond pulses of ultraviolet laser radiation, was demonstrated on ceramics [1], polymers [2,3,4,7], non-metals [5] and metals [6].

$$\Phi \left(x,y\right)=2{\Phi}_{0}\left\{cos\left[\frac{4\mathsf{\pi}x}{\lambda}sin\frac{\theta}{2}\right]+1\right\},$$

$$\Lambda =\frac{\lambda}{2sin\left(\theta /2\right)},$$

A thorough understanding of the process is essential to the precise patterning of surfaces. However, an insight into the physical mechanisms effective during laser processing cannot be gained from experimental observation, especially owing to the short laser pulse duration and the microscopic scale of the surface modification. Nevertheless, numerical simulations enable the detailed investigation of the contemplated physical effects. This approach allows for parameter variations to identify suitable process conditions with regard to texturing a specific substrate, avoiding an excessive consumption of resources in experiments. Grid-based numerical techniques, such as the finite volume and the finite element method (FEM), are usually employed in the simulation of laser material processing [9].

The mathematical model of the DLIP process originally consisted in the heat diffusion equation with the pulsed laser interference irradiation incorporated in the heat source term [8] and sink terms accommodating the latent heat of involved phase changes [10,11]. Using this model, thermal simulations of DLIP were carried out by the FEM to predict the temporal evolution of the temperature distribution near the metal surface and to assess the effect of the laser fluence on the extent of molten and vaporised material regions [10,11]. For aluminium substrates, the significant increase of the absorptivity with temperature necessitates the consideration of a temperature-dependent reflectivity in the model [12]. A detailed description of the FEM simulation of DLIP for metal substrates is presented in [13].

In this work, the thermal modelling outlined so far is expanded for the first time, according to the best of the present authors’ knowledge, to comprise the molten bath convection during nanosecond pulsed DLIP of metal surfaces. The additional complexity of this modelling approach is accepted in the prospect of insight into the role of melt pool convection in surface patterning and a profound explanation of the structuring mechanism. Furthermore, the thermofluiddynamic simulations are performed using the mesh-free smoothed particle hydrodynamics (SPH) technique. The application of mesh-free methods, which permit a deformation or even disintegration of the computational domain, is comparatively novel in the simulation of laser processes.

Specifically, the use of SPH in the modelling of laser material interactions is still little explored. The SPH method was originally developed by Gingold, Lucy and Monaghan [14,15] to address astrophysical problems. Notwithstanding its continuing importance in theoretical astrophysics, SPH was subsequently applied to various problems, notably those involving fluid flow, as in turbomachinery [16], coastal and hydraulic engineering [17]. A detailed account of these advances can be found in systematic work [18,19,20,21]. Concerning the simulation of laser processing, Chen and Beraun first solved the coupled heat conduction equations governing ultrashort laser pulse action, i.e., of subpico- to picosecond duration, on metal films by the corrective smoothed particle method [22].

The interaction of micro- to millisecond laser pulses or continuous laser irradiation with materials was modelled more frequently. Gross developed an SPH model for the laser cutting of metals [23], which was extended by Muhammad et al. to simulate the micromachining of coronary stents [24]. Considering microsecond pulses as well, Abidou et al. simulated the laser drilling of stainless steel [25]. Earlier on, Tong and Browne used weakly compressible SPH (WCSPH) to model the melt pool flow and heat transfer during laser spot welding [26], i.e., for millisecond pulses. Comparable laser spot welding simulations were performed for metallic workpieces in [27,28]. Regarding continuous irradiation, Yan et al. studied hydrodynamic interactions during laser underwater machining of alumina by SPH [29]. Later on, Hu et al. simulated conduction mode [27,30] and deep penetration laser welding [30] of aluminium. Russell et al. developed a comprehensive SPH methodology for laser based additive manufacturing and applied it to selective laser track melting [28]. Tanaka et al. investigated also heat conduction due to a stationary or moving laser source using the moving particle semi-implicit method similar to SPH [31]. In addition to [28], selective laser melting was modelled by WCSPH in [32,33,34].

On the contrary, SPH was rarely employed to address the effects of nanosecond laser pulses. The authors of this manuscript previously suggested a thermal model of DLIP for metallic substrates [35]. Cao and Shin predicted the particle motion due to phase explosion during high fluence laser ablation of metals by SPH [36]. Further, Alshaer et al. used SPH to simulate the thermal ablation of aluminium at elevated fluences, particularly the ejection of particles by the recoil pressure [37]. However, the melt pool flow during nanosecond laser irradiation at moderate fluences was not studied to date using SPH. According to the best of the authors’ knowledge, an incompressible SPH (ISPH) model was not applied to the laser-induced molten bath flow before, in contrast to weld pool convection [38].

This section presents a mathematical model of the material behaviour during single pulse DLIP. The equations governing the considered physical phenomena are stated and subsequently rewritten in non-dimensional form to reveal the dimensionless numbers characterising the process.

Throughout the DLIP process, the energy conservation is of fundamental importance. Correspondingly, the mixed enthalpy-temperature formulation of the heat transfer equation reads

$$\rho \frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}h}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}t}=\kappa \Delta T+{\dot{q}}^{\u2034}\phantom{\rule{-0.66666pt}{0ex}}.$$

The conservation of mass and momentum is particularly significant while the substrate is molten due to the effect of the laser pulse. The associated continuity and Navier–Stokes equations are given by

$$\begin{array}{cc}\hfill \nabla \xb7v& =0,\hfill \end{array}$$

$$\begin{array}{cc}\hfill \rho \frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}v}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}t}& =-\nabla p+\eta \Delta \mathit{v}+\rho \mathit{g}.\hfill \end{array}$$

As a Lagrangian method is employed, the evolution of particle positions is governed by

$$\frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}\mathit{x}}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}t}=\mathit{v}.$$

The substantial derivative $\frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}t}=\frac{\partial}{\partial t}+\mathit{v}\xb7\nabla$ is employed on the left hand side of Equations (3), (5) and (6).

In Equation (3), the specific enthalpy h consists of sensible and latent amounts
where ${L}_{\mathrm{f}}$ and ${L}_{\mathrm{v}}$ are the latent heat of fusion and vapourisation, respectively. On the right hand side of Equation (3), $\kappa $ is the thermal conductivity, T the temperature and ${\dot{q}}^{\u2034}\phantom{\rule{-0.66666pt}{0ex}}$ denotes the heat source term. Considering two-beam interference along with a Gaussian temporal shape of the laser pulse and the Beer–Lambert absorption law, the source term is given by
where $\Phi \left(x,y\right)$ is the fluence distribution according to Equation (1), $\alpha $ is the absorption coefficient of the substrate, R is its reflectivity and $\sigma ={\tau}_{\mathrm{p}}/\left(2\sqrt{2ln2}\right)$ denotes the standard deviation of the laser pulse with the duration ${\tau}_{\mathrm{p}}$ (full width at half maximum (FWHM)) at the pulse time ${t}_{\mathrm{p}}$.

$$h={h}_{\mathrm{sens}}+{h}_{\mathrm{lat}}={c}_{\mathrm{p}}\left(T-{T}_{0}\right)+{f}_{\mathrm{m}}{L}_{\mathrm{f}}+{f}_{\mathrm{v}}{L}_{\mathrm{v}},$$

$${\dot{q}}^{\u2034}\phantom{\rule{-0.66666pt}{0ex}}\left(x,y,z,t\right)=\alpha \left(1-R\right)\frac{\Phi \left(x,y\right)}{\sigma \sqrt{2\mathsf{\pi}}}exp\left(-\frac{{\left(t-{t}_{\mathrm{p}}\right)}^{2}}{2{\sigma}^{2}}+\alpha \left(z-{z}_{\mathrm{surf}}\right)\right),$$

The molten material is conceived as an incompressible fluid, as evident from the mass and momentum conservation in Equations (4) and (5). Moreover, the Oberbeck–Boussinesq approximation is applied, i.e., the density is constant except in the body force term of the momentum Equation (5), where the density varies as a function of temperature according to
with the volumetric thermal expansion coefficient $\beta $. Further, the pressure gradient term in Equation (5) considers the total pressure comprising static and dynamic components, which is represented as
where ${p}_{\mathrm{atm}}$ is a constant atmospheric reference pressure at the material surface located at ${z}_{\mathrm{surf}}$.

$$\rho \left(T\right)={\rho}_{0}\left[1-\beta \left(T-{T}_{\mathrm{l}}\right)\right]$$

$$p={p}_{\mathrm{stat}}+{p}_{\mathrm{dyn}}={\rho}_{0}g\left({z}_{\mathrm{surf}}-z\right)+{p}_{\mathrm{atm}}+{p}_{\mathrm{dyn}},$$

Employing the wavenumber $k=2\mathsf{\pi}/\lambda $ in Equation (1), inserting Equation (1) in Equation (8), and taking the result and Equations (9) and (10) into account in Equations (3) and (5), respectively, the governing equations take the form

$$\begin{array}{cc}\hfill {\rho}_{0}\frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}h}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}t}& =\kappa \Delta T+\left(1-R\right)\frac{2{\Phi}_{0}\alpha}{\sigma \sqrt{2\mathsf{\pi}}}\left\{cos\left[2kxsin\frac{\theta}{2}\right]+1\right\}exp\left(-\frac{{\left(t-{t}_{\mathrm{p}}\right)}^{2}}{2{\sigma}^{2}}+\alpha \left(z-{z}_{\mathrm{surf}}\right)\right),\hfill \\ \hfill \nabla \xb7\mathit{v}& =0,\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\rho}_{0}\frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}\mathit{v}}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}t}& =-\nabla {p}_{\mathrm{dyn}}+\eta \Delta \mathit{v}-\beta \left(T-{T}_{\mathrm{l}}\right){\rho}_{0}\mathit{g},\hfill \\ \hfill \frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}\mathit{x}}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}t}& =\mathit{v}.\hfill \end{array}$$

Due to the short time scale of nanosecond laser interference patterning, heat losses due to convection and radiation are neglected and the heat transfer equation is subject to adiabatic boundary conditions

$$-\kappa \frac{\partial T}{\partial \mathit{n}}=0,\phantom{\rule{1.em}{0ex}}\frac{\partial h}{\partial \mathit{n}}=0.$$

In the presence of melt, homogeneous Neumann boundary conditions are applied to the dynamic pressure and the no-slip condition, i.e., vanishing velocity, is enforced at the bottom of the molten pool

$$\frac{\partial {p}_{\mathrm{dyn}}}{\partial \mathit{n}}=0,\phantom{\rule{1.em}{0ex}}\mathit{v}=0.$$

At the free surface, an inhomogeneous Neumann condition, the Marangoni boundary condition, is employed for the horizontal velocity, cf. [39,40], and a vanishing vertical velocity is prescribed

$$\eta \frac{\partial u}{\partial z}=\frac{\partial \gamma}{\partial x}=\frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}\gamma}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}T}\frac{\partial T}{\partial x},\phantom{\rule{1.em}{0ex}}w=0.$$

To rewrite the governing equations in dimensionless form, the non-dimensionalisation of the variables position, time, velocity, pressure, temperature and specific enthalpy is performed employing the scales L, ${L}^{2}/a$, $a/L$, ${\rho}_{0}{a}^{2}/{L}^{2}$, ${T}_{\mathrm{v}}-{T}_{0}$ and ${c}_{\mathrm{p},\mathrm{ref}}\left({T}_{\mathrm{v}}-{T}_{0}\right)$, respectively. In particular, the characteristic length scale is specified as the diffusion length $L=2\sqrt{a{\tau}_{\mathrm{p}}}$, where the pulse width ${\tau}_{\mathrm{p}}$ (FWHM) is chosen as the laser beam dwell time. Applying the aforementioned scales, the resulting dimensionless variables are

$${\mathit{x}}^{*}=\frac{\mathit{x}}{2\sqrt{a{\tau}_{\mathrm{p}}}},\phantom{\rule{4pt}{0ex}}{t}^{*}=\frac{t}{4{\tau}_{\mathrm{p}}},\phantom{\rule{4pt}{0ex}}{\mathit{v}}^{*}=2\sqrt{\frac{{\tau}_{\mathrm{p}}}{a}}\mathit{v},\phantom{\rule{4pt}{0ex}}{p}_{\mathrm{dyn}}^{*}=\frac{4{\tau}_{\mathrm{p}}}{{\rho}_{0}a}{p}_{\mathrm{dyn}},\phantom{\rule{4pt}{0ex}}{T}^{*}=\frac{T-{T}_{0}}{{T}_{\mathrm{v}}-{T}_{0}},\phantom{\rule{4pt}{0ex}}{h}^{*}=\frac{h}{{c}_{\mathrm{p},\mathrm{ref}}\left({T}_{\mathrm{v}}-{T}_{0}\right)}.$$

Using the dimensionless variables in Equation (16), Equations (4), (6), (11) and (12) are obtained in dimensionless form as
with the dimensionless standard deviation ${\sigma}^{*}=\sigma /\left(4{\tau}_{\mathrm{p}}\right)$ of the laser pulse, periodicity ${\Lambda}^{*}=\Lambda /L$ and absorption coefficient ${\alpha}^{*}=\alpha L$ in Equation (17). In particular, the dimensionless form of the Marangoni boundary condition in Equation (15) is given as [40]

$$\begin{array}{cc}\hfill \frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{h}^{*}}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{t}^{*}}& =\Delta {T}^{*}+\left(1-R\right)\frac{\mathit{La}}{{\sigma}^{*}\sqrt{2\mathsf{\pi}}}\left\{cos\left[\frac{2\mathsf{\pi}{x}^{*}}{{\Lambda}^{*}}\right]+1\right\}exp\left(-\frac{{\left({t}^{*}-{t}_{\mathrm{p}}^{*}\right)}^{2}}{2{\sigma}^{*\phantom{\rule{0.166667em}{0ex}}2}}+{\alpha}^{*}\left({z}^{*}-{z}_{\mathrm{surf}}^{*}\right)\right),\hfill \end{array}$$

$$\begin{array}{cc}\hfill \nabla \xb7{\mathit{v}}^{*}& =0,\hfill \end{array}$$

$$\begin{array}{cc}\hfill \frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{\mathit{v}}^{*}}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{t}^{*}}& =-\nabla {p}_{\mathrm{dyn}}^{*}+\mathit{Pr}\Delta {\mathit{v}}^{*}+\mathit{Pr}\mathit{Ra}\frac{{T}^{*}-{T}_{\mathrm{l}}^{*}}{1-{T}_{\mathrm{l}}^{*}}{\mathit{e}}_{z},\hfill \end{array}$$

$$\begin{array}{cc}\hfill \frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{\mathit{x}}^{*}}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{t}^{*}}& ={\mathit{v}}^{*}\hfill \end{array}$$

$$\frac{\partial {u}^{*}}{\partial {z}^{*}}=-\frac{\mathit{Ma}}{1-{T}_{\mathrm{l}}^{*}}\frac{\partial {T}^{*}}{\partial {x}^{*}}.$$

The dimensionless numbers emerging in the dimensionless Equations (17)–(21), complemented by the Fourier number $\mathit{Fo}$ corresponding to the considered physical time, are the Laser number $\mathit{La}$
the Prandtl number $\mathit{Pr}$, the Rayleigh number $\mathit{Ra}$ and the Marangoni number $\mathit{Ma}$ defined by

$$\mathit{La}=\frac{2{\Phi}_{0}\alpha}{{\rho}_{\mathrm{ref}}{c}_{\mathrm{p},\mathrm{ref}}\left({T}_{\mathrm{v}}-{T}_{0}\right)},\phantom{\rule{1.em}{0ex}}\mathit{Fo}=\frac{{t}_{\mathrm{end}}}{4{\tau}_{\mathrm{p}}},$$

$$\mathit{Pr}=\frac{\nu}{a},\phantom{\rule{1.em}{0ex}}\mathit{Ra}=8{\tau}_{\mathrm{p}}\frac{\beta \left({T}_{\mathrm{v}}-{T}_{\mathrm{l}}\right)g\sqrt{a{\tau}_{\mathrm{p}}}}{\nu},\phantom{\rule{1.em}{0ex}}\mathit{Ma}=-\frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}\gamma}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}T}\frac{{T}_{\mathrm{v}}-{T}_{\mathrm{l}}}{{\rho}_{0}\nu}2\sqrt{\frac{{\tau}_{\mathrm{p}}}{a}}.$$

In addition, the dimensionless specific enthalpy is obtained as
where the phase change numbers of melting and vapourisation are defined as

$${h}^{*}={T}^{*}+{f}_{\mathrm{m}}{\mathit{Ph}}_{\mathrm{s}/\mathrm{l}}+{f}_{\mathrm{v}}{\mathit{Ph}}_{\mathrm{l}/\mathrm{v}},$$

$${\mathit{Ph}}_{\mathrm{s}/\mathrm{l}}=\frac{{L}_{\mathrm{f}}}{{c}_{\mathrm{p},\mathrm{ref}}\left({T}_{\mathrm{v}}-{T}_{0}\right)},\phantom{\rule{1.em}{0ex}}{\mathit{Ph}}_{\mathrm{l}/\mathrm{v}}=\frac{{L}_{\mathrm{v}}}{{c}_{\mathrm{p},\mathrm{ref}}\left({T}_{\mathrm{v}}-{T}_{0}\right)}.$$

Finally, the molten and vaporised mass fractions required for Equation (24) and the relation between enthalpy and temperature are given by

$$\begin{array}{cc}\hfill {f}_{\mathrm{m}}& =\left\{\begin{array}{cc}0& {h}^{*}\le {T}_{\mathrm{s}}^{*}\hfill \\ \frac{{h}^{*}-{T}_{\mathrm{s}}^{*}}{{T}_{\mathrm{l}}^{*}-{T}_{\mathrm{s}}^{*}+{\mathit{Ph}}_{\mathrm{s}/\mathrm{l}}}& {T}_{\mathrm{s}}^{*}<{h}^{*}<{T}_{\mathrm{l}}^{*}+{\mathit{Ph}}_{\mathrm{s}/\mathrm{l}}\hfill \\ 1& {T}_{\mathrm{l}}^{*}+{\mathit{Ph}}_{\mathrm{s}/\mathrm{l}}\le {h}^{*}\hfill \end{array},\right.\hfill \end{array}$$

$$\begin{array}{cc}\hfill {f}_{\mathrm{v}}& =\left\{\begin{array}{cc}0& {h}^{*}\le 1+{\mathit{Ph}}_{\mathrm{s}/\mathrm{l}}\hfill \\ \frac{{h}^{*}-1-{\mathit{Ph}}_{\mathrm{s}/\mathrm{l}}}{{\mathit{Ph}}_{\mathrm{l}/\mathrm{v}}}& 0<{h}^{*}-1-{\mathit{Ph}}_{\mathrm{s}/\mathrm{l}}<{\mathit{Ph}}_{\mathrm{l}/\mathrm{v}}\hfill \\ 1& 1+{\mathit{Ph}}_{\mathrm{s}/\mathrm{l}}+{\mathit{Ph}}_{\mathrm{l}/\mathrm{v}}\le {h}^{*}\hfill \end{array},\right.\hfill \end{array}$$

$$\begin{array}{cc}\hfill {T}^{*}& =\left\{\begin{array}{cc}{h}^{*}& {f}_{\mathrm{m}}={f}_{\mathrm{v}}=0\hfill \\ {T}_{\mathrm{s}}^{*}+{f}_{\mathrm{m}}\left({T}_{\mathrm{l}}^{*}-{T}_{\mathrm{s}}^{*}\right)& 0<{f}_{\mathrm{m}}<1,{f}_{\mathrm{v}}=0\hfill \\ {h}^{*}-{\mathit{Ph}}_{\mathrm{s}/\mathrm{l}}& {f}_{\mathrm{m}}=1,{f}_{\mathrm{v}}=0\hfill \\ {T}_{\mathrm{v}}^{*}=1& {f}_{\mathrm{m}}=1,0<{f}_{\mathrm{v}}\le 1\hfill \end{array}.\right.\hfill \end{array}$$

The mesh-free SPH method employed in the present work is described in the following. The exposition of the method starts with its fundamentals and then explains the important aspects to be considered to set up a working SPH algorithm. This account comprises the incompressible SPH (ISPH) approach pursued to solve the fluid flow along with a discussion of numerical stability.

The mathematical identity underlying the SPH method is the representation of a quantity $\phi $ by its convolution with the delta distribution $\delta $ given as

$$\phi \left(\mathit{x}\right)=\int \phi \left({\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right)\delta \left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right)\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}.$$

The kernel approximation is performed by replacing the delta distribution inside the integral in Equation (29) with a kernel function W, resulting in
where ${\Omega}_{\mathit{x}}$ denotes the compact support of W centred at $\mathit{x}$ and defined by the smoothing length ${l}_{\mathrm{sm}}$. The integration in Equation (30) is carried out as a summation over discrete particles representing the computational domain $\Omega $, giving rise to the particle approximation

$$\phi \left(\mathit{x}\right)\approx {\int}_{{\Omega}_{\mathit{x}}}\phi \left({\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right)W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},$$

$$\phi \left(\mathit{x}\right)\approx \sum _{j=1}^{N}\frac{{m}_{j}}{{\rho}_{j}}\phi \left({\mathit{x}}_{j}\right)W\left(\mathit{x}-{\mathit{x}}_{j},{l}_{\mathrm{sm}}\right).$$

The requirements for the kernel function imposed to establish the approximation in Equation (30) and further desirable properties are [41,42]

$$\begin{array}{ccc}\hfill \underset{{l}_{\mathrm{sm}}\to 0}{lim}W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)& =\delta \left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right)\hfill & \hfill \left(\mathrm{approximation}\text{}\mathrm{of}\text{}\delta \text{}\mathrm{distribution}\right),\end{array}$$

$$\begin{array}{ccc}\hfill W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)& \ge 0\hfill & \hfill \text{\hspace{1em}}\forall {\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\in {\Omega}_{\mathit{x}}\phantom{\rule{2.em}{0ex}}\phantom{\rule{2.em}{0ex}}\phantom{\rule{2.em}{0ex}}\phantom{\rule{2.em}{0ex}}\phantom{\rule{2.em}{0ex}}\phantom{\rule{0.277778em}{0ex}}\left(\mathrm{positivity}\right),\end{array}$$

$$\begin{array}{ccc}\hfill W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)& =0\hfill & \hfill \text{\hspace{1em}}\forall {\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\notin {\Omega}_{\mathit{x}}\phantom{\rule{2.em}{0ex}}\phantom{\rule{2.em}{0ex}}\phantom{\rule{2.em}{0ex}}\phantom{\rule{1.em}{0ex}}\left(\mathrm{compact}\text{}\mathrm{support}\right),\end{array}$$

$$\begin{array}{ccc}\hfill W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)& =W\left(\left|\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right|,{l}_{\mathrm{sm}}\right)\hfill & \hfill \left(\mathrm{spherical}\text{}\mathrm{symmetry}\right),\end{array}$$

$$\begin{array}{ccc}\hfill W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)& \ge W\left(\mathit{x}-{\mathit{x}}^{\u2033}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)\hfill & \hfill \forall \left|\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right|<\left|\mathit{x}-{\mathit{x}}^{\u2033}\phantom{\rule{-0.66666pt}{0ex}}\right|\phantom{\rule{2.em}{0ex}}\phantom{\rule{1.em}{0ex}}\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}\left(\mathrm{monotonicity}\right),\end{array}$$

$$\begin{array}{ccc}\hfill {\int}_{{\Omega}_{\mathit{x}}}W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}& =1\hfill & \hfill \left(\mathrm{normalisation}\right),\end{array}$$

$$\begin{array}{ccc}\hfill W& \in {C}^{2}\hfill & \hfill \left(\mathrm{smoothness}\right).\end{array}$$

The existence of a continuous second derivative of the kernel function in condition (38) is imposed [43], and, in particular, is necessary if it is evaluated in second derivative approximations.

In the present work, a quintic B-spline ascribed to Schoenberg [44], introduced in SPH context by Morris et al. [45], is employed. In the notation used by Speith [46], this kernel function consisting of piecewise quintic polynomials is given as
where $r=\left|\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right|$, the problem dimension is denoted by d and the normalisation constant results from Equation (37) as [46]

$$W\left(r,{l}_{\mathrm{sm}}\right)=\frac{\varsigma}{{l}_{\mathrm{sm}}^{d}}\left\{\begin{array}{cc}{\left(1-r/{l}_{\mathrm{sm}}\right)}^{5}-6{\left(\frac{2}{3}-r/{l}_{\mathrm{sm}}\right)}^{5}+15{\left(\frac{1}{3}-r/{l}_{\mathrm{sm}}\right)}^{5}& 0\le r/{l}_{\mathrm{sm}}<\frac{1}{3},\hfill \\ {\left(1-r/{l}_{\mathrm{sm}}\right)}^{5}-6{\left(\frac{2}{3}-r/{l}_{\mathrm{sm}}\right)}^{5}& \frac{1}{3}\le r/{l}_{\mathrm{sm}}<\frac{2}{3},\hfill \\ {\left(1-r/{l}_{\mathrm{sm}}\right)}^{5}& \frac{2}{3}\le r/{l}_{\mathrm{sm}}<1,\hfill \\ 0& 1\le r/{l}_{\mathrm{sm}},\hfill \end{array}\right.$$

$$\varsigma =\left\{\begin{array}{cc}243/40& d=1,\hfill \\ 15309/\left(478\mathsf{\pi}\right)& d=2,\hfill \\ 2187/\left(40\mathsf{\pi}\right)& d=3.\hfill \end{array}\right.$$

Unlike the representation of the physical quantity $\phi $ itself in Equation (31), the particle approximation of a derivative term involves the first derivative of the kernel function, as elucidated in Section 3.2. In case of the quintic spline kernel function presented above in Equation (39), the first derivative is given as

$$\frac{\partial W\left(r,{l}_{\mathrm{sm}}\right)}{\partial r}=\frac{-5\varsigma}{{l}_{\mathrm{sm}}^{d+1}}\left\{\begin{array}{cc}{\left(1-r/{l}_{\mathrm{sm}}\right)}^{4}-6{\left(\frac{2}{3}-r/{l}_{\mathrm{sm}}\right)}^{4}+15{\left(\frac{1}{3}-r/{l}_{\mathrm{sm}}\right)}^{4}& 0\le r/{l}_{\mathrm{sm}}<\frac{1}{3},\hfill \\ {\left(1-r/{l}_{\mathrm{sm}}\right)}^{4}-6{\left(\frac{2}{3}-r/{l}_{\mathrm{sm}}\right)}^{4}& \frac{1}{3}\le r/{l}_{\mathrm{sm}}<\frac{2}{3},\hfill \\ {\left(1-r/{l}_{\mathrm{sm}}\right)}^{4}& \frac{2}{3}\le r/{l}_{\mathrm{sm}}<1,\hfill \\ 0& 1\le r/{l}_{\mathrm{sm}}.\hfill \end{array}\right.$$

An approximation of the gradient of a scalar field $\phi $ can be obtained by inserting $\nabla \phi $ in Equation (30) and performing an integration by parts [19]

$$\begin{array}{cc}\hfill \nabla \phi \left(\mathit{x}\right)& \approx {\int}_{\Omega \cap {\Omega}_{\mathit{x}}}{\nabla}_{\phantom{\rule{-1.99997pt}{0ex}}{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}}\left[\phi \left({\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right)\right]W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& ={\int}_{\Omega \cap {\Omega}_{\mathit{x}}}{\nabla}_{\phantom{\rule{-1.99997pt}{0ex}}{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}}\left[\phi \left({\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right)W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)\right]\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}-{\int}_{\Omega \cap {\Omega}_{\mathit{x}}}\phi \left({\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right){\nabla}_{\phantom{\rule{-1.99997pt}{0ex}}{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\phantom{\rule{0.166667em}{0ex}}}W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& ={\oint}_{\partial \Omega \cap {\Omega}_{\mathit{x}}}\phi \left({\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right)W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)\mathit{n}\left({\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right)\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}\Gamma +{\int}_{\Omega \cap {\Omega}_{\mathit{x}}}\phi \left({\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right){\nabla}_{\phantom{\rule{-1.99997pt}{0ex}}\mathit{x}\phantom{\rule{0.166667em}{0ex}}}W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}.\hfill \end{array}$$

Concerning the right hand side of Equation (42), the first integral is rewritten as a surface integral according to Gauss’ theorem. In addition, the symmetry of the kernel function in Equation (35) implies the property ${\nabla}_{\phantom{\rule{-1.99997pt}{0ex}}{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\phantom{\rule{0.166667em}{0ex}}}W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)=-{\nabla}_{\phantom{\rule{-1.99997pt}{0ex}}\mathit{x}\phantom{\rule{0.166667em}{0ex}}}W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)$ of the kernel gradient, which is applied to the second integral.

As the surface integral in Equation (42) vanishes if $\mathit{x}$ is far enough from the domain boundary $\partial \Omega $, the straightforward approximation of the gradient of a scalar quantity $\phi $ results as

$$\nabla \phi \left(\mathit{x}\right)\approx {\int}_{{\Omega}_{\mathit{x}}}\phi \left({\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\right){\nabla}_{\phantom{\rule{-1.99997pt}{0ex}}\mathit{x}\phantom{\rule{0.166667em}{0ex}}}W\left(\mathit{x}-{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}},{l}_{\mathrm{sm}}\right)\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}{\mathit{x}}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}.$$

Consequently, the discrete particle approximation of the gradient of a scalar field $\phi $ is given as

$$\nabla \phi \left(\mathit{x}\right)\approx \sum _{j=1}^{N}\frac{{m}_{j}}{{\rho}_{j}}\phi \left({\mathit{x}}_{j}\right)\frac{\mathit{x}-{\mathit{x}}_{j}}{\left|\mathit{x}-{\mathit{x}}_{j}\right|}\frac{\partial W\left(\mathit{x}-{\mathit{x}}_{j},{l}_{\mathrm{sm}}\right)}{\partial \left|\mathit{x}-{\mathit{x}}_{j}\right|}.$$

It is remarkable that the derivative of a quantity, e.g., the gradient of a scalar field $\phi $ in Equation (44), is evaluated using the quantity at the particle positions and the known derivative of the kernel function, see Equation (41), in the SPH method. However, if the discrete gradient presented in Equation (44) is applied to a constant field, the resulting value is not equal to zero [19].

Intending to solve a system of partial differential equations using a mesh-free particle method, the discrete differential operators are evaluated at the position ${\mathit{x}}_{i}$ of a particle of interest i in the following. The zeroth order consistency of the gradient approximation can be recovered by employing the symmetric gradient operator [19,47]
where the shorthand notation ${\phi}_{i}=\phi \left({\mathit{x}}_{i}\right)$ is used. On the contrary, it is not recommended to approximate the pressure gradient in the momentum Equation (5) by the symmetric gradient operator in Equation (45), as ${\mathit{G}}_{i}^{-}$ does not fulfil the action–reaction principle [47,48], i.e., Newton’s third law. Linear momentum conservation can be ensured using the antisymmetric gradient operator [42,48]
which satisfies the action–reaction principle.

$${\mathit{G}}_{i}^{-}\left({\phi}_{j}\right)=\frac{1}{{\rho}_{i}}\sum _{j=1}^{{N}_{i}}{m}_{j}\left({\phi}_{j}-{\phi}_{i}\right)\frac{{\mathit{x}}_{i}-{\mathit{x}}_{j}}{\left|{\mathit{x}}_{i}-{\mathit{x}}_{j}\right|}\frac{\partial W\left({\mathit{x}}_{i}-{\mathit{x}}_{j},{l}_{\mathrm{sm}}\right)}{\partial \left|{\mathit{x}}_{i}-{\mathit{x}}_{j}\right|}\approx {\left(\nabla \phi \right)}_{i},$$

$${\mathit{G}}_{i}^{+}\left({\phi}_{j}\right)={\rho}_{i}\sum _{j=1}^{{N}_{i}}{m}_{j}\left(\frac{{\phi}_{i}}{{\rho}_{i}^{2}}+\frac{{\phi}_{j}}{{\rho}_{j}^{2}}\right)\frac{{\mathit{x}}_{i}-{\mathit{x}}_{j}}{\left|{\mathit{x}}_{i}-{\mathit{x}}_{j}\right|}\frac{\partial W\left({\mathit{x}}_{i}-{\mathit{x}}_{j},{l}_{\mathrm{sm}}\right)}{\partial \left|{\mathit{x}}_{i}-{\mathit{x}}_{j}\right|}\approx {\left(\nabla \phi \right)}_{i},$$

In view of the accuracy of the projection-based incompressible SPH approach to be presented in Section 3.4, it is essential that the discrete gradient and divergence operators be skew-adjoint [19,47]. For this reason, the velocity divergence arising in the incompressible SPH scheme is approximated by the symmetric divergence operator [48]
where it can be shown that the discrete operators ${\mathit{G}}_{i}^{+}$ and ${D}_{i}^{-}$ are skew-adjoint [19].

$${D}_{i}^{-}\left({\mathit{v}}_{j}\right)=\frac{1}{{\rho}_{i}}\sum _{j=1}^{{N}_{i}}{m}_{j}\frac{\left({\mathit{v}}_{j}-{\mathit{v}}_{i}\right)\xb7\left({\mathit{x}}_{i}-{\mathit{x}}_{j}\right)}{\left|{\mathit{x}}_{i}-{\mathit{x}}_{j}\right|}\frac{\partial W\left({\mathit{x}}_{i}-{\mathit{x}}_{j},{l}_{\mathrm{sm}}\right)}{\partial \left|{\mathit{x}}_{i}-{\mathit{x}}_{j}\right|}\approx {\left(\nabla \xb7\mathit{v}\right)}_{i},$$

Moreover, a second-order differential operator, the Laplacian, is required to solve the governing equations, e.g., for the approximation of the heat conduction term in the heat transfer Equation (3) and the viscous diffusion term in the momentum Equation (5). If a procedure analogous to the one in Equation (42) is followed, the resulting expression involves the second derivative of the kernel function [19,46,47]. Therefore, the approximation is very sensitive to particle disorder [48] and plagued by an undetermined sign of the summands [19,49] as the kernel function exhibits a point of inflexion. On the other hand, a discrete Laplacian could be constructed by composing a gradient and a divergence operator [19,50], i.e., ${\Delta}_{i}\approx {D}_{i}^{-}{\mathit{G}}_{j}^{+}$. This exact operator comprises a double summation, which makes it impracticable due to the high computational effort [19,47].

Nevertheless, an approximate Laplacian can be obtained, as first suggested by Morris et al. [45], by combining a finite difference expression for the gradient with an SPH divergence operator [47]. The resulting second-order differential operator employs only the first derivative of the kernel function. This discrete Laplacian originates in the modelling of heat conduction [51,52] and viscous diffusion [45]. Considering the approximation of a general diffusion term, the Laplacian is given by [45,52]
where ${\Gamma}_{\phi}$ denotes the diffusion coefficient related to the quantity $\phi $. The Laplacian ${L}_{i}$ in Equation (48) represents a simplified form for a spherically symmetric kernel function. Instead of a scalar $\phi $, the approximate Laplacian can also be applied to a vectorial quantity and the resulting vector is denoted by ${\mathit{L}}_{i}$ in this case.

$${L}_{i}\left({\Gamma}_{\phi ,j},{\phi}_{j}\right)=\sum _{j=1}^{{N}_{i}}\frac{{m}_{j}}{{\rho}_{j}}\left({\Gamma}_{\phi ,i}+{\Gamma}_{\phi ,j}\right)\frac{{\phi}_{i}-{\phi}_{j}}{\left|{\mathit{x}}_{i}-{\mathit{x}}_{j}\right|}\frac{\partial W\left({\mathit{x}}_{i}-{\mathit{x}}_{j},{l}_{\mathrm{sm}}\right)}{\partial \left|{\mathit{x}}_{i}-{\mathit{x}}_{j}\right|}\approx {\left\{\nabla \xb7\left({\Gamma}_{\phi}\nabla \phi \right)\right\}}_{i},$$

The treatment of boundary conditions in SPH received little attention from its inception in theoretical astrophysics, where boundaries do not play a crucial role [53]. Nevertheless, the application of boundary conditions in an SPH algorithm is a nontrivial task. Particles close to the boundary exhibit an incomplete kernel support due to the truncation by the boundary, the particle deficiency problem [54].

Different approaches were proposed to remedy the incompleteness of the kernel support near the boundaries, including three classical boundary treatment techniques. The first is the generation of mirror particles by direct reflexion of near-wall particles across the boundary, as proposed by Libersky and Petschek [55]. On the other hand, Monaghan [56] introduced boundary particles, which are located at the edges of the computational domain and exert repulsive forces on approaching particles. In addition, Randles and Libersky [54] completed the kernel support of near-wall particles by fixed dummy particles situated beyond the boundaries.

Furthermore, Kulasegaram et al. [57] devised a profound treatment for wall boundaries using a variational formulation of the SPH equations of motion. In this approach, a correction factor, i.e., a function characterising the completeness of the kernel support, for particles near a wall is introduced and, on its basis, boundary contact forces are evaluated [57]. Considering the distance of a particle from the nearest boundary segment, the correction factor and its derivative are determined using a spline approximation of the original kernel integral [19,57]. Different methods were also developed to account for the case of an intersection of boundary segments [57].

In the present work, boundary conditions to the energy Equation (3) are imposed using fixed dummy particles located beyond the boundaries. Moreover, a combination of edge particles situated right on the boundary and dummy particles is employed to apply the boundary conditions required for the solution of molten pool convection using the incompressible SPH method presented in the following Section 3.4, as suggested by different authors [58,59].

The traditional approach of treating incompressible fluid motion using SPH consists of the assumption of a slight compressibility of the fluid and the numerical solution of the compressible conservation equations. In addition, this weakly compressible SPH (WCSPH) method comprises an equation of state to close the system of the compressible continuity and Navier–Stokes equations. However, this procedure gives rise to considerable noise in the pressure field, since small fluctuations in the density field result in large pressure fluctuations due to the stiffness of the equation of state [59]. Remedies suggested for this disadvantage of WCSPH include a particle initialisation algorithm [60] and the introduction of diffusive corrections [61,62] in combination with a particle shifting approach [63].

More recently than WCSPH, the projection method for the numerical solution of the incompressible Navier–Stokes equations developed independently by Chorin [64] and Temam [65] was introduced in SPH context by Cummins and Rudman [50]. The projection method employs the Helmholtz–Hodge decomposition theorem which states that any vector field $\mathit{w}$ on a smoothly bounded domain $\Omega $ can be uniquely decomposed in the form
where $\mathit{u}$ is a solenoidal vector field and parallel to the boundary $\partial \Omega $, and $\nabla p$ is an irrotational vector field [66].

$$\mathit{w}=\mathit{u}+\nabla p,$$

The SPH approach relying on the projection method for solving the equations governing incompressible fluid flow belongs to the incompressible SPH (ISPH) methods. Meanwhile, several ISPH algorithms based on the projection method were proposed using different formulations of the incompressibility constraint, i.e., imposing either a divergence-free velocity field [50] or density invariance [58] or combining both criteria [67]. Other than the projection-based approach, an ISPH method using Lagrange multipliers acting as non-thermodynamic pressures to enforce constant particle volume was presented in [68].

It was recognised that the sole requirement of zero velocity divergence in a projection-based ISPH algorithm leads to substantial particle density variations [67], due to the stretching and compression of particle positions [69]. This particle clustering phenomenon arises due to the ordered motion of the particles along the streamlines [19,69]. By enforcing density invariance in the projection-based ISPH algorithm, the problem can be avoided [67], but this method exhibits a reduced numerical accuracy [69]. Combining both zero velocity divergence and density invariance requirements in an ISPH approach, Hu and Adams could remedy the excessive particle density variations [67].

The latter method being associated with an increased computational effort; Xu et al. proposed an alternative [69]. This approach combines the projection-based ISPH method enforcing a divergence-free velocity field with a slight shift of the particle positions across the streamlines [69]. As this ISPH algorithm is expected to produce acceptable results with reasonable computational effort, it is employed in the present work. For the numerical solution of the system of governing Equations (4), (6), (11) and (12) in the time step ${t}^{n+1}={t}^{n}+\Delta t$, the steps of the algorithm are given below. (The gradient of a vector field in Equation (57) is a second-order tensor equal to the Jacobian, i.e., $\nabla {\mathit{v}}_{i}^{n+1}={\left[\nabla {v}_{1,i}^{n+1}\cdots \nabla {v}_{d,i}^{n+1}\right]}^{\mathsf{T}}$.)

$$\begin{array}{ccc}\hfill \overline{\mathit{x}}& ={\mathit{x}}^{n}+{\mathit{v}}^{n}\Delta t\hfill & \hfill \text{\hspace{1em}}\mathrm{position}\text{}\mathrm{advection}\end{array}$$

$$\begin{array}{ccc}\hfill \overline{\mathit{v}}& ={\mathit{v}}^{n}+\left(\nu \Delta {\mathit{v}}^{n}+\beta \left({T}^{n}-{T}_{\mathrm{l}}\right)g{\mathit{e}}_{z}\right)\Delta t\hfill & \hfill \text{\hspace{1em}}\mathrm{velocity}\text{}\mathrm{prediction}\end{array}$$

$$\begin{array}{ccc}\hfill \Delta {p}_{\mathrm{dyn}}^{n+1}& =\frac{{\rho}_{0}}{\Delta t}\nabla \xb7\overline{\mathit{v}}\hfill & \hfill \text{\hspace{1em}}\mathrm{solution}\text{}\mathrm{of}\text{}\mathrm{pressure}\text{}\mathrm{Poisson}\text{}\mathrm{equation}\text{}\left(\mathrm{PPE}\right)\end{array}$$

$$\begin{array}{ccc}\hfill {\mathit{v}}^{n+1}& =\overline{\mathit{v}}-\frac{\Delta t}{{\rho}_{0}}\nabla {p}_{\mathrm{dyn}}^{n+1}\hfill & \hfill \text{\hspace{1em}}\mathrm{velocity}\text{}\mathrm{correction}\end{array}$$

$$\begin{array}{ccc}\hfill {\mathit{x}}^{n+1}& ={\mathit{x}}^{n}+\left({\mathit{v}}^{n}+{\mathit{v}}^{n+1}\right)\frac{\Delta t}{2}\hfill & \hfill \text{\hspace{1em}}\mathrm{position}\text{}\mathrm{update}\end{array}$$

$$\begin{array}{ccc}\hfill {h}^{n+1}& ={h}^{n}+\left(\kappa \Delta {T}^{n+1}+{\dot{q}}^{\u2034}{\phantom{\rule{-0.66666pt}{0ex}}}^{,n+1}\right)\frac{\Delta t}{{\rho}_{0}}\hfill & \hfill \text{\hspace{1em}}\mathrm{enthalpy}\text{}\mathrm{update}\end{array}$$

$$\begin{array}{ccc}\hfill {\tilde{\mathit{x}}}_{i}^{\phantom{\rule{0.166667em}{0ex}}n+1}\phantom{\rule{-0.166667em}{0ex}}& ={\mathit{x}}_{i}^{n+1}+C\alpha {\mathit{R}}_{i}\hfill & \hfill \text{\hspace{1em}}\mathrm{position}\text{}\mathrm{shift}\end{array}$$

$$\begin{array}{ccc}\hfill {\tilde{\mathit{v}}}_{i}^{\phantom{\rule{0.166667em}{0ex}}n+1}\phantom{\rule{-0.166667em}{0ex}}& ={\mathit{v}}_{i}^{n+1}+\nabla {\mathit{v}}_{i}^{n+1}\left({\tilde{\mathit{x}}}_{i}^{\phantom{\rule{0.166667em}{0ex}}n+1}-{\mathit{x}}_{i}^{n+1}\right)\hfill & \hfill \text{\hspace{1em}}\mathrm{velocity}\text{}\mathrm{adjustment}\end{array}$$

$$\begin{array}{ccc}\hfill {\tilde{h}}_{i}^{n+1}& ={h}_{i}^{n+1}+{\left(\nabla {h}_{i}^{n+1}\right)}^{\mathsf{T}}\left({\tilde{\mathit{x}}}_{i}^{\phantom{\rule{0.166667em}{0ex}}n+1}-{\mathit{x}}_{i}^{n+1}\right)\hfill & \hfill \text{\hspace{1em}}\mathrm{enthalpy}\text{}\mathrm{adjustment}\end{array}$$

The time step size has to satisfy several constraints in order to ensure the numerical stability of the SPH algorithm. In short, the time step can be determined according to
the Courant–Friedrichs–Lewy, maximum particle acceleration, viscous and thermal diffusion, and surface tension conditions on the time step being given within brackets.

$$\Delta t=min\left(0.25\frac{{l}_{\mathrm{sm}}}{{v}_{max}},0.25\underset{i}{min}\sqrt{\frac{{l}_{\mathrm{sm}}}{\left|{\mathit{f}}_{i}\right|}},0.125\underset{i}{min}\frac{{l}_{\mathrm{sm}}^{2}}{{\nu}_{i}},0.125\underset{i}{min}\frac{{l}_{\mathrm{sm}}^{2}}{{a}_{i}},0.25\sqrt{\frac{\rho {l}_{\mathrm{sm}}^{3}}{2\mathsf{\pi}\gamma}}\right),$$

To identify the interacting pairs of neighbouring particles, the cell index method [70] is employed. This approach was described by Hockney and Eastwood [71] with regard to the evaluation of short-range forces in particle methods and first used in the SPH method by Monaghan and Gingold [72]. The principal idea is to subdivide the domain into square cells with a side length equal to or slightly larger than the interaction radius [70,71]. Each particle is assigned to a cell on the basis of its position. The cells are represented by an array of linked lists maintained to keep track of the particles residing in each cell [70,71,72]. Therefore, the search for all neighbours of a given particle reduces to the examination of the particles in the same cell and in the eight adjacent cells (in 2D) [70,71]. The number of tests required for determining all interacting particle pairs is further halved in the present case of symmetric interactions. The computational effort may be further decreased using an optimal, larger cell size [73] or combining the cell linked lists with a Verlet list, which introduces additional memory requirements [74].

The intention of the present work is to perform numerical simulations of heat transfer and fluid flow during single pulse DLIP using SPH. A 2D section of the substrate in the $x-z$ plane comprising the interaction zone due to one period of the interference pattern is considered for this purpose, see Figure 1. This area is discretised using particles as illustrated in Section 4.1. Thereafter, the concrete numerical scheme employed to solve the dimensionless governing Equations (17)–(20) is presented.

The basic strategy to solve the energy Equation (17) in the computational domain throughout the simulation duration is shown in Section 4.2. Shortly after the onset of the nanosecond laser pulse, a thin layer of material adjacent to the surface starts to melt in the vicinity of the interference maximum. In the subdomain representing the molten pool, the solution of the energy equation is a part of a more comprehensive approach to solve the complete set of Equations (17)–(20), which is clarified in Section 4.3.

As indicated above, the present research considers a rectangular computational domain in the $x-z$ plane, which is discretised using particles. The length of the domain is given by the period $\Lambda $ of the interference pattern in Equation (2), whereas its height amounts to several diffusion lengths L, as defined in Section 2.2. To maintain a reasonable computational effort, the employed particle distribution exhibits a local refinement towards the surface in line with the strategy the authors used before in [35]. In particular, the interaction zone, where the substrate is expected to melt due to the action of the laser pulse, is discretised by equidistant fine particles. Despite the associated high computational effort, this idea is followed to allow for equal size particles in the subdomain representing the molten pool.

In detail, the discretisation is performed starting from the bottom of the computational domain, where a row of coarse particles of 1 µm diameter is employed. The discretisation is continued using a successive reduction of the particle size towards the surface, the diameter of the particles in the rows situated above being given by a geometric sequence with a common ratio, or quotient, of $q=7/9$. A graphical representation of the discretisation is provided in Figure 2, where Figure 2a shows the coarser part referred to so far. As mentioned in the foregoing paragraph, the refinement does not go beyond a minimum particle size, and numerous rows of particles of this size are arranged on a Cartesian lattice to discretise the interaction zone adjacent to the surface, see Figure 2b. The uniform and initially equidistant distribution of fine particles is employed to avoid potential detrimental effects of different particle sizes during the numerical solution of molten pool convection. As a trade-off between an appropriate resolution of the absorption length and a manageable computational effort, the minimum particle diameter is specified as ${\left(7/9\right)}^{22}$ µm $\approx 3.97$ nm.

In addition, the discretisation is extended by three rows and columns of dummy particles situated beyond the horizontal and vertical domain boundaries, respectively. Consequently, the kernel function support of the adjacent interior particles is completed by the dummy particles, with the respective particle diameter being defined by the one of the nearest interior particle. This provision of dummy particles is prescribed by the employed smoothing length. The latter is related to the (dimensionless) kernel support radius ${r}_{W}$, which coincides with the number of segments of the radial coordinate in the definition of the spline kernel function in Equation (39).

As stated by Morris, the number of interacting particles should be augmented for a kernel function with larger compact support [75]. For the commonly employed cubic spline kernel function, the kernel support radius is ${r}_{W}=2$ and a typical smoothing length amounts to ${l}_{\mathrm{sm}}^{*}=2.4\Delta {x}^{*}$. As the quintic spline kernel function exhibits a larger compact support with ${r}_{W}=3$, an extended smoothing length
is used here, where $\Delta {x}^{*}$ is the initial separation of particles on a Cartesian grid. The smoothing length specified in Equation (60) corresponds to a neighbourhood of 45 interacting particles in an equidistant rectangular arrangement with particle spacing $\Delta {x}^{*}$ in $d=2$ dimensions.

$${l}_{\mathrm{sm}}^{*}=3.75\Delta {x}^{*}$$

Note that for interactions between particles of different size, the smoothing length is averaged as explained in the following. Consider two interacting particles i and j separated by $k\in \mathbb{N}\phantom{\rule{-0.166667em}{0ex}}\setminus \phantom{\rule{-0.166667em}{0ex}}\left\{0\right\}$ refinements given by a geometric sequence with common ratio $q<1$, and assume without loss of generality that the diameter ${d}_{i}^{*}$ of particle i is larger than the diameter ${d}_{j}^{*}={d}_{i}^{*}{q}^{k}$ of particle j. The vertical distance between these two particles can be written as
i.e., the vertical separation between these particles is the k-fold average vertical particle distance

$${z}_{j}^{*}-{z}_{i}^{*}=k\Delta {z}_{ji}^{*,\mathrm{av}}=\frac{{d}_{i}^{*}}{2}\left(1+q\right)\sum _{l=0}^{k-1}{q}^{l}=\frac{{d}_{i}^{*}}{2}\left(1+q\right)\frac{1-{q}^{k}}{1-q},$$

$$\Delta {z}_{ji}^{*,\mathrm{av}}=\frac{{d}_{i}^{*}}{2k}\left(1+q\right)\frac{1-{q}^{k}}{1-q}=\frac{{d}_{j}^{*}}{2k}\left(1+{q}^{-1}\right)\frac{1-{q}^{-k}}{1-{q}^{-1}}.$$

From Equation (60), it is evident that the strength of interaction between particles i and j is non-vanishing only for $k\le 3$. In particular, the average vertical particle spacing in Equation (62) reduces to the arithmetic mean of the particle diameters in the case $k=1$. The above consideration leads to the averaged smoothing length used for the interaction between particles i and j
where the smoothing length ${l}_{\mathrm{sm},i}^{*}$ for particle i is given by Equation (60) with the particle separation being equal to the particle diameter, i.e., ${\left(\Delta {x}^{*}\right)}_{i}={d}_{i}^{*}$.

$${l}_{\mathrm{sm},ij}^{*}=\frac{{l}_{\mathrm{sm},i}^{*}}{2k}\left(1+q\right)\frac{1-{q}^{k}}{1-q}=\frac{{l}_{\mathrm{sm},i}^{*}}{2k}\left(1+q\right)\sum _{l=0}^{k-1}{q}^{l}={l}_{\mathrm{sm},ji}^{*},$$

The energy Equation (17) is solved using the methodology presented earlier by the authors and their co-authors [35]. In particular, the energy Equation (17) in mixed enthalpy–temperature formulation is implicitly integrated in time. Consequently, the discretisation of the dimensionless energy Equation (17) in time for an interior or top edge particle i is written as
where the approximate Laplacian from Equation (48) and the power of the laser heat source per unit mass ${\dot{Q}}_{i}^{*,n+1}$ are employed in dimensionless form.

$$\frac{{h}_{i}^{*,n+1}-{h}_{i}^{*,n}}{\Delta {t}^{*}}={L}_{i}^{*}\left(1,{T}_{j}^{*,n+1}\right)+{\dot{Q}}_{i}^{*,n+1},$$

In the discrete Laplacian ${L}_{i}^{*}(1,{T}_{j}^{*,n+1})$, the temperature ${T}_{i}^{*,n+1}$ inside the summation is rewritten as a function of the specific enthalpy ${h}_{i}^{*,n+1}$ according to Equation (24). Subsequently, Equation (64) can be rearranged to result in an expression for the specific enthalpy at the new time step [35]

$${h}_{i}^{*,n+1}=\frac{{h}_{i}^{*,n}+\Delta {t}^{*}\left({\dot{Q}}_{i}^{*,n+1}-2{\sum}_{j=1}^{{N}_{i}}\frac{{m}_{j}^{*}}{{\rho}_{j}^{*}}\frac{{T}_{j}^{*,n+1}+{f}_{\mathrm{m},i}^{n+1}{\mathit{Ph}}_{\mathrm{s}/\mathrm{l}}+{f}_{\mathrm{v},i}^{n+1}{\mathit{Ph}}_{\mathrm{l}/\mathrm{v}}}{{r}_{ij}^{*,n+1}}\frac{\partial W\left({r}_{ij}^{*,n+1},{l}_{\mathrm{sm},ij}^{*}\right)}{\partial {r}_{ij}^{*,n+1}}\right)}{1-2\Delta {t}^{*}{\sum}_{j=1}^{{N}_{i}}\frac{{m}_{j}^{*}}{{\rho}_{j}^{*}}\frac{1}{{r}_{ij}^{*,n+1}}\frac{\partial W\left({r}_{ij}^{*,n+1},{l}_{\mathrm{sm},ij}^{*}\right)}{\partial {r}_{ij}^{*,n+1}}}.$$

Replacing also ${T}_{j}^{*,n+1}$ in Equation (65) with the equivalent terms from Equation (24), a linear system of equations is obtained for the enthalpy field at the new time step, which is constituted by the individual equations for all interior and top edge particles i given as [35]

$$\begin{array}{cc}\hfill \phantom{\rule{1.em}{0ex}}& \left(1-2\Delta {t}^{*}\sum _{j=1}^{{N}_{i}}\frac{{m}_{j}^{*}}{{\rho}_{j}^{*}}\frac{1}{{r}_{ij}^{*,n+1}}\frac{\partial W\left({r}_{ij}^{*,n+1},{l}_{\mathrm{sm},ij}^{*}\right)}{\partial {r}_{ij}^{*,n+1}}\right){h}_{i}^{*,n+1}+2\Delta {t}^{*}\sum _{j=1}^{{N}_{i}}\frac{{m}_{j}^{*}}{{\rho}_{j}^{*}}\frac{{h}_{j}^{*,n+1}}{{r}_{ij}^{*,n+1}}\frac{\partial W\left({r}_{ij}^{*,n+1},{l}_{\mathrm{sm},ij}^{*}\right)}{\partial {r}_{ij}^{*,n+1}}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& +2\Delta {t}^{*}\sum _{j=1}^{{N}_{i}}\frac{{m}_{j}^{*}}{{\rho}_{j}^{*}}\frac{\left({f}_{\mathrm{m},i}^{n+1}-{f}_{\mathrm{m},j}^{n+1}\right){\mathit{Ph}}_{\mathrm{s}/\mathrm{l}}+\left({f}_{\mathrm{v},i}^{n+1}-{f}_{\mathrm{v},j}^{n+1}\right){\mathit{Ph}}_{\mathrm{l}/\mathrm{v}}}{{r}_{ij}^{*,n+1}}\frac{\partial W\left({r}_{ij}^{*,n+1},{l}_{\mathrm{sm},ij}^{*}\right)}{\partial {r}_{ij}^{*,n+1}}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& \phantom{\rule{2.em}{0ex}}={h}_{i}^{*,n}+\Delta {t}^{*}\left(1-R\right)\frac{\mathit{La}}{{\sigma}^{*}\sqrt{2\mathsf{\pi}}}\left\{cos\left[\frac{2\mathsf{\pi}{x}_{i}^{*,n+1}}{{\Lambda}^{*}}\right]+1\right\}{e}^{-\frac{{\left({t}^{*,n+1}-{t}_{\mathrm{p}}^{*}\right)}^{2}}{2{\sigma}^{*\phantom{\rule{0.166667em}{0ex}}2}}+{\alpha}^{*}\left({z}_{i}^{*,n+1}-{z}_{\mathrm{surf}}^{*}\right)}.\hfill \end{array}$$

Concerning the iterative solution of this linear system of equations, a few aspects of interest are given in the following. To begin with, the third term on the left hand side of Equation (66) is considered not to contribute to the system matrix, i.e., the molten and vaporised mass fractions are not conceived as a function of the unknown enthalpy in the present iteration. Instead, this term is calculated based on the determination of the molten and vaporised mass fractions from the previous enthalpy iterate.

The formulation in Equation (66) gives rise to a linear system of equations characterised in the following. The system matrix is symmetric, given particles of equal volume (or equal mass in a uniform density approach), its main diagonal elements are positive and the matrix is strictly diagonally dominant. These three properties imply that the matrix is positive definite [76]. The conjugate gradient (CG) method introduced by Hestenes and Stiefel [77] is commonly used to iteratively solve a linear system with a symmetric and positive definite matrix.

In the present work, a preconditioned variant of the CG algorithm [78] in conjunction with a Jacobi preconditioner is employed for the iterative solution of the linear system arising from the discretised dimensionless heat transfer Equation (64). However, only a single CG step is performed, then the quantities in Equations (26)–(28) are updated according to the new enthalpy iterate and this iterative procedure is restarted.

As indicated in the introductory paragraph of this section, the subsection at hand provides a detailed account of the ISPH scheme employed to solve the system of dimensionless Equations (17)–(20) in a manner analogous to the algorithm given in Equations (50)–(58) for the numerical treatment of Equations (4), (6), (11) and (12). Due to the intricacy of this approach, the subsection is further subdivided for the sake of clarity. The numerical details of the individual steps, notably the respective discrete differential operators employed, of the ISPH algorithm are presented in Section 4.3.1. Aspects to be considered for the numerical solution of the pressure Poisson equation (PPE), which plays a pivotal role in the projection-based ISPH method, are covered in Section 4.3.2.

While the substrate is locally molten as a result of the thermalisation of the interference irradiation provided by the laser pulse, the melt flow is calculated using the ISPH algorithm explained below. The particle positions and velocities are evolved for the completely fluid particles inside the molten pool, whereas the dynamic pressure values are also determined for the surrounding edge particles. In accordance with the sole solution of the energy Equation (17) presented in Section 4.2, the transition from the old time step ${t}^{*,n}$ to the new time step ${t}^{*,n+1}={t}^{*,n}+\Delta {t}^{*}$ is considered for the numerical treatment of Equations (17)–(20) shown here.

At first, the fluid particles are advected to intermediate positions based on the old velocity

$${\overline{\mathit{x}}}_{i}^{*}={\mathit{x}}_{i}^{*,n}+{\mathit{v}}_{i}^{*,n}\Delta {t}^{*}.$$

Considering the acceleration due to viscous and body forces, an intermediate velocity field
is predicted, which is not divergence-free. Subsequently, the PPE for enforcing zero velocity divergence
is solved at the fluid particles inside and the edge particles surrounding the molten pool. The intermediate velocity is then corrected using the gradient of the determined dynamic pressure field to obtain a divergence-free velocity field

$${\overline{\mathit{v}}}_{i}^{*}={\mathit{v}}_{i}^{*,n}+\left(\mathit{Pr}{\mathit{L}}_{i}^{*}\left(1,{\mathit{v}}_{j}^{*,n}\right)+\mathit{Pr}\mathit{Ra}\frac{{T}_{i}^{*,n}-{T}_{\mathrm{l}}^{*}}{1-{T}_{\mathrm{l}}^{*}}{\mathit{e}}_{z}\right)\Delta {t}^{*}$$

$${L}_{i}^{*}\left(1,{p}_{\mathrm{dyn},j}^{*,n+1}\right)=\frac{1}{\Delta {t}^{*}}{D}_{i}^{-\phantom{\rule{0.166667em}{0ex}}*}\phantom{\rule{-0.166667em}{0ex}}\left({\overline{\mathit{v}}}_{j}^{*}\right)$$

$${\mathit{v}}_{i}^{*,n+1}={\overline{\mathit{v}}}_{i}^{*}-\Delta {t}^{*}{\mathit{G}}_{i}^{+\phantom{\rule{0.166667em}{0ex}}*}\phantom{\rule{-0.166667em}{0ex}}\left({p}_{\mathrm{dyn},j}^{*,n+1}\right).$$

After that, the particle positions are updated using both old and new velocity fields

$${\mathit{x}}_{i}^{*,n+1}={\mathit{x}}_{i}^{*,n}+\left({\mathit{v}}_{i}^{*,n}+{\mathit{v}}_{i}^{*,n+1}\right)\frac{\Delta {t}^{*}}{2}.$$

The specific enthalpy at the new time step is calculated
as discussed above in Section 4.2. To avoid a too orderly particle motion along the streamlines, the particle positions are slightly shifted [69]
with a constant $C\in \left[0.01,0.1\right]$, the shifting magnitude ${\alpha}^{*}={max}_{j}\left|{\mathit{v}}_{j}^{*,n+1}\right|\Delta {t}^{*}$ and the shifting vector ${\mathit{R}}_{i}$ depending on the average particle spacing ${r}_{i}^{*,\mathrm{av}}$ and the unit vectors ${\mathit{n}}_{ij}$. A truncated Taylor series expansion is then employed to adjust the velocity values to the final particle positions (The discrete gradient operator in Equation (74) is characterised by ${\mathit{G}}_{i}^{-\phantom{\rule{0.166667em}{0ex}}*}\left({\mathit{v}}_{j}^{*,n+1}\right)={\left[{\mathit{G}}_{i}^{-\phantom{\rule{0.166667em}{0ex}}*}\phantom{\rule{-0.166667em}{0ex}}\left({v}_{x,j}^{*,n+1}\right)\phantom{\rule{4pt}{0ex}}{\mathit{G}}_{i}^{-\phantom{\rule{0.166667em}{0ex}}*}\phantom{\rule{-0.166667em}{0ex}}\left({v}_{z,j}^{*,n+1}\right)\right]}^{\mathsf{T}}$.)

$${h}_{i}^{*,n+1}={h}_{i}^{*,n}+\Delta {t}^{*}\left({L}_{i}^{*}\left(1,{T}_{j}^{*,n+1}\right)+{\dot{Q}}_{i}^{*,n+1}\right)$$

$${\tilde{\mathit{x}}}_{i}^{\phantom{\rule{0.166667em}{0ex}}*,n+1}\phantom{\rule{-0.166667em}{0ex}}={\mathit{x}}_{i}^{*,n+1}+C{\alpha}^{*}{\mathit{R}}_{i},\phantom{\rule{1.em}{0ex}}\phantom{\rule{4pt}{0ex}}{\mathit{R}}_{i}=\sum _{j=1}^{{N}_{i}}{\left(\frac{{r}_{i}^{*,\mathrm{av}}}{{r}_{ij}^{*,n+1}}\right)}^{2}\phantom{\rule{-0.166667em}{0ex}}{\mathit{n}}_{ij},\phantom{\rule{1.em}{0ex}}\phantom{\rule{4pt}{0ex}}{r}_{i}^{*,\mathrm{av}}=\frac{1}{{N}_{i}}\sum _{j=1}^{{N}_{i}}{r}_{ij}^{*,n+1},\phantom{\rule{1.em}{0ex}}\phantom{\rule{4pt}{0ex}}{\mathit{n}}_{ij}=\frac{{\mathit{x}}_{ij}^{*,n+1}}{{r}_{ij}^{*,n+1}},$$

$${\tilde{\mathit{v}}}_{i}^{\phantom{\rule{0.166667em}{0ex}}*,n+1}\phantom{\rule{-0.166667em}{0ex}}={\mathit{v}}_{i}^{*,n+1}+{\mathit{G}}_{i}^{-\phantom{\rule{0.166667em}{0ex}}*}\phantom{\rule{-0.166667em}{0ex}}\left({\mathit{v}}_{j}^{*,n+1}\right)\left({\tilde{\mathit{x}}}_{i}^{\phantom{\rule{0.166667em}{0ex}}*,n+1}\phantom{\rule{-0.166667em}{0ex}}-{\mathit{x}}_{i}^{*,n+1}\right).$$

In addition, an analogous adjustment is performed for the specific enthalpy values

$${\tilde{{h}_{i}}}^{*,n+1}={h}_{i}^{*,n+1}+{\left({\mathit{G}}_{i}^{-\phantom{\rule{0.166667em}{0ex}}*}\phantom{\rule{-0.166667em}{0ex}}\left({h}_{j}^{*,n+1}\right)\right)}^{\mathsf{T}}\left({\tilde{\mathit{x}}}_{i}^{\phantom{\rule{0.166667em}{0ex}}*,n+1}\phantom{\rule{-0.166667em}{0ex}}-{\mathit{x}}_{i}^{*,n+1}\right).$$

Finally, the adjusted specific enthalpy values are considered for the determination of the molten and vaporised mass fractions and the particle temperatures ${T}_{i}^{*,n+1}$ according to Equations (26)–(28).

It is observed that the procedure described above uses two different sets of positions for the particle approximations during each time step. The discrete differential operators employed in the velocity prediction, solution of the PPE and velocity correction steps in Equations (68)–(70) rely on the advected particle positions given in Equation (67). On the other hand, the discrete summations in the specific enthalpy update, position correction, velocity and specific enthalpy adjustment steps in Equations (72)–(75) are based on the updated particle positions calculated in Equation (71).

In addition, after the specific enthalpy and temperature update in Equation (72), the temperature gradient is determined using the symmetric gradient operator ${\mathit{G}}_{i}^{-\phantom{\rule{0.166667em}{0ex}}*}\phantom{\rule{-0.166667em}{0ex}}\left({T}_{j}^{*,n+1}\right)$, the horizontal component being required for the Marangoni boundary condition given by Equation (21) in the next time step.

Furthermore, the ISPH scheme exposed above is subject to restrictions on the dimensionless time step size. Rewriting Equation (59) in dimensionless form, the conditions to be respected are given as
where the Ohnesorge number $\mathit{Oh}=\nu /\sqrt{\gamma L/\rho}$ emerges in the dimensionless surface tension condition.

$$\Delta {t}^{*}=min\left(0.25\frac{{l}_{\mathrm{sm}}^{*}}{{v}_{max}^{*}},0.25\underset{i}{min}\sqrt{\frac{{l}_{\mathrm{sm}}^{*}}{\left|{\mathit{f}}_{i}^{*}\right|}},0.125\underset{i}{min}\frac{{l}_{\mathrm{sm}}^{*\phantom{\rule{0.166667em}{0ex}}2}}{\mathit{Pr}},0.125\underset{i}{min}{l}_{\mathrm{sm}}^{*\phantom{\rule{0.166667em}{0ex}}2},0.25\frac{\mathit{Oh}}{\mathit{Pr}}\sqrt{\frac{{l}_{\mathrm{sm}}^{*\phantom{\rule{0.166667em}{0ex}}3}}{2\mathsf{\pi}}}\right),$$

Applying the dimensionless form of the discrete Laplacian and symmetric divergence operator in Equations (47) and (48), respectively, the zero velocity divergence PPE (69) to be solved is written out as

$$2\sum _{j=1}^{{N}_{i}}\frac{{m}_{j}^{*}}{{\rho}_{j}^{*}}\frac{{p}_{\mathrm{dyn},i}^{*,n+1}-{p}_{\mathrm{dyn},j}^{*,n+1}}{\left|{\overline{\mathit{x}}}_{i}^{*}-{\overline{\mathit{x}}}_{j}^{*}\right|}\frac{\partial W\left({\overline{\mathit{x}}}_{i}^{*}-{\overline{\mathit{x}}}_{j}^{*},{l}_{\mathrm{sm}}^{*}\right)}{\partial \left|{\overline{\mathit{x}}}_{i}^{*}-{\overline{\mathit{x}}}_{j}^{*}\right|}=\frac{1}{\Delta {t}^{*}}\sum _{j=1}^{{N}_{i}}\frac{{m}_{j}^{*}}{{\rho}_{j}^{*}}\frac{\left({\overline{\mathit{v}}}_{j}^{*}-{\overline{\mathit{v}}}_{i}^{*}\right)\xb7\left({\overline{\mathit{x}}}_{i}^{*}-{\overline{\mathit{x}}}_{j}^{*}\right)}{\left|{\overline{\mathit{x}}}_{i}^{*}-{\overline{\mathit{x}}}_{j}^{*}\right|}\frac{\partial W\left({\overline{\mathit{x}}}_{i}^{*}-{\overline{\mathit{x}}}_{j}^{*},{l}_{\mathrm{sm}}^{*}\right)}{\partial \left|{\overline{\mathit{x}}}_{i}^{*}-{\overline{\mathit{x}}}_{j}^{*}\right|}.$$

The presence of Equation (77), which all fluid particles inside and edge particles bounding the molten pool are to satisfy, constitutes a linear system of equations $\mathit{A}{\mathit{p}}_{\mathrm{dyn}}^{*,n+1}=\mathit{b}$ for the dynamic pressure values.

It should be evident from the explanations given below that the concept of a system matrix $\mathit{A}$, although it is not assembled in the numerical code, is meaningful for the solution of the PPE. The dynamic pressure field to be determined by solving Equation (77) is subject to a homogeneous Neumann boundary condition, the dimensionless form of the first condition in Equation (14), at the molten pool edges, as mentioned in Section 2.1. For this type of boundary condition, the technique adopted in this work requires that the dynamic pressure values at the edge particles are assigned to the respective dummy particles situated beyond the molten pool edges in the outward normal direction.

The idea is to write the linear system consisting of an equation of the type given in Equation (77) for each fluid and edge particle i in the form $\mathit{A}{\mathit{p}}_{\mathrm{dyn}}^{*,n+1}=\mathit{b}$ with a square matrix $\mathit{A}$, where ${\mathit{p}}_{\mathrm{dyn}}^{*,n+1}$ is the solution vector comprising the dynamic pressures of all fluid and edge particles i. Contemplating on the left hand side of Equation (77) for fluid and edge particles, the structure of the system matrix can be understood from the following observations. To begin with, as the derivative of the kernel function vanishes at the origin, the self-interaction of a fluid or edge particle i is neglected.

The notation $\left\{{i}_{\mathrm{D}}\right\}$ is introduced for the set of dummy particles associated with an edge particle i. In particular, the interactions of an edge particle i with the dummy particles $j\in \left\{{i}_{\mathrm{D}}\right\}$ are disregarded as the dynamic pressure values cancel out. Apart from the interaction with a (different) edge particle j, each contribution due to the interaction of a fluid or edge particle i with any dummy particle ${j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\in \left\{{j}_{\mathrm{D}}\right\}$ by means of the value ${p}_{\mathrm{dyn},j}^{*,n+1}$ is allocated to the associated edge particle j with equal dynamic pressure. These considerations lead to the non-zero entries of the matrix $\mathit{A}$ listed in Table 1. In detail, the off-diagonal elements ${a}_{ij}$ in the second row of Table 1 are valid for both fluid and edge particles i.

The following properties of the system matrix $\mathit{A}$ can be inferred from its elements presented in Table 1 and further reflexions. In general, the matrix is comparatively sparse since the number of interactions ${N}_{i}$ of neighbouring particles with particle i is much less than the total number of particles representing the molten pool and its edges. Nevertheless, the situation that each fluid particle can in principle interact with any other fluid or edge particle throughout the simulation leads to the fact that the matrix is not banded. It is evident from Table 1 that the main diagonal entries ${a}_{ii}$ of the matrix $\mathit{A}$ are negative, whereas its off-diagonal elements ${a}_{ij}$ are positive.

The matrix $\mathit{A}$ is weakly diagonally dominant, as $\forall i:\left|{a}_{ii}\right|={\sum}_{j\ne i}\left|{a}_{ij}\right|$. More precisely, the sum of the elements in each row of $\mathit{A}$ is zero. Therefore, the vector ${\mathit{x}}_{0}={\left(1,1,\dots ,1\right)}^{\mathsf{T}}$ is in the kernel of $\mathit{A}$. Consequently, the matrix $\mathit{A}$ is not of full rank, i.e., $\mathit{A}$ is singular, according to the rank-nullity theorem. This problematic feature of $\mathit{A}$ arises for confined flows in the absence of a Dirichlet boundary condition on the pressure [47,79].

Now consider a fluid particle i situated sufficiently close to an edge so that there is at least one dummy particle ${j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}$ beyond edge particle j interacting with particle i, i.e., $\exists {j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\in \left\{{j}_{\mathrm{D}}\right\}:W\left({\overline{r}}_{i{j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}}^{*},{l}_{\mathrm{sm}}^{*}\right)>0$, then ${a}_{ij}=-2\frac{{m}_{j}^{*}}{{\rho}_{j}^{*}}\frac{1}{{\overline{r}}_{ij}^{*}}\frac{\partial W\left({\overline{r}}_{ij}^{*},{l}_{\mathrm{sm}}^{*}\right)}{\partial {\overline{r}}_{ij}^{*}}-2{\sum}_{{j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\in \left\{{j}_{\mathrm{D}}\right\}}\frac{{m}_{{j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}}^{*}}{{\rho}_{{j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}}^{*}}\frac{1}{{\overline{r}}_{i{j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}}^{*}}\frac{\partial W\left({\overline{r}}_{i{j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}}^{*},{l}_{\mathrm{sm}}^{*}\right)}{\partial {\overline{r}}_{i{j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}}^{*}}$. On the contrary, there is no additional interaction allocated to the one with fluid particle i in the row giving rise to Equation (77) for edge particle j, i.e., ${a}_{ji}=-2\frac{{m}_{i}^{*}}{{\rho}_{i}^{*}}\frac{1}{{\overline{r}}_{ji}^{*}}\frac{\partial W\left({\overline{r}}_{ji}^{*},{l}_{\mathrm{sm}}^{*}\right)}{\partial {\overline{r}}_{ji}^{*}}$. Consequently, as ${a}_{ij}\ne {a}_{ji}$, the matrix $\mathit{A}$ is non-symmetric.

As the system matrix $\mathit{A}$ is not regular as mentioned above, it should be regularised to solve the linear system of equations. This is performed through a slight reinforcement of the diagonal entries and by subtracting from the right hand side of the equation system the mean of its entries [47].

In this work, the linear system with a resulting non-symmetric and regular system matrix is solved iteratively using the BiCGStab algorithm [80]. However, a numerical instability of BiCGStab denoted as pivot breakdown [81] is often encountered during the solution of the PPE. For this reason, a stabilisation of the algorithm is implemented, which detects a pivot near-breakdown situation and consequently restarts the BiCGStab iteration.

In the following, the model presented above is applied to investigate DLIP of metallic substrates using a single nanosecond pulse. For the sake of clarity, the present section is divided into two subsections. First, the parameters of the DLIP process are indicated in Section 5.1 along with the material properties considered for both stainless steel and aluminium substrates as well as the actual values of the dimensionless numbers. Subsequently, the details of the numerical investigation and the results of DLIP simulations by means of (I)SPH are presented in Section 5.2.

As already indicated above, the interference of two coherent laser beams giving rise to a sinusoidal intensity distribution is studied in this work. The considered parameters of the DLIP process, particularly with regard to the laser heat source, are given in Table 2. Note that the periodicity $\Lambda $ is given as a function of the laser wavelength $\lambda $ and the angle of intersection $\theta $ between the interfering beams according to Equation (2). Due to the specification of the thermal diffusion length as the characteristic length scale, the Fourier number defined in Equation (22) and given in Table 2 depends only on the pulse duration and the considered physical time, irrespective of the substrate.

Concerning the considered substrates, AISI 304 stainless steel and high-purity aluminium, the relevant material properties are listed in Table 3. If available, the respective temperature-dependent values are averaged over the entire interval from the room temperature to the vapourisation point to obtain the reference density ${\rho}_{\mathrm{ref}}$, specific heat ${c}_{\mathrm{p},\mathrm{ref}}$, thermal conductivity ${\kappa}_{\mathrm{ref}}$ and thermal diffusivity ${a}_{\mathrm{ref}}$. Based on the latter, the thermal diffusion length L is calculated and, along with process parameters and material properties from Table 2 and Table 3, employed to determine the dimensionless quantities stated in Table 4.

Furthermore, in case of the stainless steel substrate, the absence of surfactants, e.g., sulphur, in the alloy is assumed and a constant temperature coefficient of surface tension $\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}\gamma /\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}T$ is considered, see Table 3. The nonlinear temperature dependency of both surface tension and its temperature coefficient [82] in the presence of surface-active elements as well as the effect of this additional complexity on DLIP results of stainless steel substrates are explored in a separate investigation.

In particular, the absorptivity of aluminium exhibits a very distinct temperature dependency. Although the reflectivity of aluminium is 92.5% at the wavelength $\lambda =$ 355 nm at room temperature [85,95] as indicated in Table 3, it decreases to 65% at the melting point [84]. Therefore, the effect of temperature on the reflectivity of aluminium is considered by employing the following relation [12]

$$R=\left\{\begin{array}{cc}1.048-3.989\xb7{10}^{-4}T/\mathrm{K}& {T}_{0}\le T<{T}_{\mathrm{m}},\hfill \\ 0.65& {T}_{\mathrm{m}}\le T.\hfill \end{array}\right.$$

Simulations of DLIP experiments using a single nanosecond pulse are performed by means of the approach presented above. The fluence of the interference pattern to be specified allows local melting of the substrate at the interference maxima without significant vapourisation of material. To this end, the suitable range of fluences is identified from the literature [96,97] and this processing parameter is reported in Table 5. In addition, the Laser number and the depth of the near-surface zone, which is discretised by the finest particles initially arranged on an equidistant Cartesian lattice, provided for the molten pool computation, are indicated in Table 5 for the different fluences and metallic substrates.

A few comments on the principal characteristics of the DLIP simulations are given beforehand. Due to the absorption and thermalisation of the incident radiation, the substrate is rapidly heated after the onset of the laser pulse and locally melted. The surface temperature at the interference maximum is saved in each time step of the SPH simulation to understand its evolution. In the case of an elevated laser fluence and a high absorptivity, the maximum surface temperature may even reach the vapourisation point. As opposed to the previous work [35], the vapourisation of particles by their respective deactivation is disregarded in the present investigation, although the simulations suggest that the latent heat of vapourisation is completely inserted into surface particles for elevated energy densities starting from 0.6 J/cm^{2} for stainless steel and from 1.0 J/cm^{2} for aluminium.

In the presence of molten metal, the nonuniform surface temperature gives rise to surface tension variations, the highest values being found at the boundaries of the molten pool in accordance with the negative temperature coefficient of surface tension. The surface tension gradients cause shear stresses at the surface inducing an outward flow from the centre of the melt pool, i.e., at the interference maximum, to its edge and, owing to continuity, a backward flow at its bottom. To get an idea of the importance of melt pool convection, the melt pool dimensions and the extremes of the horizontal particle velocity found below the surface particles are recorded in each time step of the SPH simulation while metal melt is present. After the action of the laser pulse, the temperature field inside the substrate tends to homogenise as a result of heat conduction. This phenomenon in particular entails the resolidification of the molten material, starting from the melt pool edges.

Considering the stainless steel substrate, SPH simulations of single pulse DLIP are performed for laser energy densities ranging from 0.4 J/cm^{2} to 1.0 J/cm^{2}. The predictions of the time-dependent surface temperature at the interference maximum are presented in Figure 3 for this range of fluences along with the temporal intensity variation of the laser pulse. Except for the lowest fluence, the surface is heated to the vapourisation temperature at the interference maximum as indicated in the figure.

The melt pool width and depth obtained from the SPH computations are presented in Figure 4a,b, respectively, for DLIP of stainless steel in the considered range of energy densities. For plotting the trends in Figure 4, each discrete value of a melt pool dimension is employed only once at an intermediate point in time to avoid a stair-step appearance of the graphs, particularly in Figure 4b. It is observed from Figure 4 that the presence of the melt pool is considerably prolonged for an enhanced laser fluence. In addition, the melt pool width and even more its depth increase distinctly with the energy density for low and moderate fluences up to 0.7 J/cm^{2}. Nevertheless, the melt pool developing due to the effect of the laser pulse is rather shallow with an aspect ratio (depth/width) of the order of 0.1.

The resulting maximum (outward) horizontal velocity trends during melt pool convection are reported in Figure 5 for the investigated laser fluences. The velocity variations in Figure 5 show a superlinear growth of the maximum horizontal velocity just below the surface for an increased laser energy density. In particular, the velocity trends exhibit an excessive increase of the maximum horizontal velocity magnitude starting at a time later than 55 ns, i.e., after the laser pulse duration (FWHM), except for the simulation considering the lowest fluence.

The latter observation is examined in detail for a fluence of 0.5 J/cm^{2} in Figure 6, where the evolution of the temperature, the horizontal temperature gradient and the velocity at the melt pool surface is presented in the relevant temporal range. After the onset of melt cooling, temperature gradients as high as several 1000 K/µm, see Figure 6b, arise at particles in the close proximity of a central superficial zone at the vapourisation point, see Figure 6a. By means of the Marangoni boundary condition in Equation (15), these excessive temperature gradients lead to a pair of inner horizontal velocity extrema reaching a magnitude of 4 m/s, in addition to the pre-existing pair of outer horizontal velocity extrema with absolute values not much higher than 3.5 m/s, see Figure 6c. As the inner velocity maximum temporarily exceeds the outer one, it is the former that is reported in Figure 5 during this time interval. However, for the moderate fluence considered, the inner velocity extrema eventually decay in the absence of a surface zone at the vapourisation point and the adjacent high temperature gradients.

In addition, to provide some more detail to the horizontal velocity trends shown in Figure 5, the temporal evolution of the magnitude of inner and outer horizontal velocity maxima is presented in Figure 7 for different laser fluences. For a moderate fluence of 0.5 J/cm^{2}, Figure 7a shows that the outer and inner extrema coexist for more than 4 ns, which is also observed in Figure 6c. However, the magnitude of the inner velocity extrema exceeds the one at the outer extrema only for 1.7 ns according to Figure 7a, which is evident from Figure 5 as well. On the other hand, the horizontal velocity magnitudes presented in Figure 7b reveal that the inner extrema emerge within fractions of a nanosecond and soon dominate for elevated laser energy densities. In this situation, the (outer) horizontal velocity extrema existing from the onset of molten pool convection subsequently vanish, cf. Figure 7b.

The results of DLIP simulations performed by SPH for the high-purity aluminium substrate and energy densities ranging from 0.6 J/cm^{2} to 1.2 J/cm^{2} are outlined in the following. At first, Figure 8 presents the temporal evolution of the maximum surface temperature for different fluences. Moderate temperatures of the aluminium surface at the interference maximum are reported in Figure 8, even at elevated fluences. This result can be attributed to the low absorptivity, particularly in the solid state, the high thermal diffusivity and the enhanced specific heat and latent heat of fusion of aluminium. The trends in Figure 8 further indicate a rapid homogenisation of the substrate temperature after the laser pulse action, the maximum surface temperature returning to the melting point comparatively early.

The melt pool dimensions, i.e., its width and depth, predicted by the SPH simulations are plotted against time in Figure 9a,b, respectively. For the graphs in Figure 9, the individual discrete values of the melt pool dimensions are again considered only once at a central point in time to evade a stair-step presentation. The figure indicates that an increase in the laser fluence results in an extended presence as well as larger dimensions of the melt pool. Concerning DLIP of aluminium, the data presented in Figure 9b reveal that the melt pool is comparatively deep, where an aspect ratio of the melt pool in the range of 0.2 to 0.3 can be deduced from Figure 9.

Furthermore, the maximum horizontal velocity magnitudes at particles below the surface are presented in Figure 10 for the considered laser fluences. The velocity trends displayed in Figure 10 confirm a superlinear growth of the horizontal velocity below the surface for an increased laser energy density. In fact, remarkable velocity magnitudes beyond 10 m/s are attained even at moderate laser fluences employed for DLIP of aluminium substrates. This observation is explored in greater detail in Figure 11, where the evolution of the temperature, the horizontal temperature gradient and the velocity are presented along the melt pool surface and the adjacent solid material.

More precisely, the aforementioned quantities are displayed in Figure 11 for DLIP of aluminium using an energy density of 1.0 J/cm^{2} in the temporal range characterised by the onset of cooling and the passing of the maximum melt pool width and horizontal velocity, cf. Figure 9a and Figure 10. At the beginning of the temperature field homogenisation, Figure 11b shows that high temperature gradients arise in the close vicinity of a central surface zone at the vapourisation point, cf. Figure 11a. However, these thermal gradients do not result in the complete formation of separate horizontal velocity extrema adjacent to the central high temperature region, see Figure 11c. Apart from or in the absence of this surface region, Figure 11b exposes a consistent temperature gradient on both sides of the interference maximum. Owing to the temperature dependence of surface tension, this steady temperature gradient along the melt pool surface gives rise to a smoothly varying horizontal velocity profile with its extrema apart from the centre and closer to the edges of the melt pool, see Figure 11c.

The results presented above suggest that stainless steel and aluminium represent two very dissimilar substrate materials with distinct behaviour when subject to laser interference irradiation. To begin with, the maximum surface temperature exhibits a more pronounced sensitivity to the laser fluence for the aluminium substrate, see Figure 8, as compared to the results in Figure 3 for stainless steel. This observation may be attributed to the high reflectivity of aluminium and the self-enhancing effect of its temperature-dependent absorptivity during laser heating. In addition, the thermal diffusivity of aluminium, being one order of magnitude higher than the one of stainless steel, facilitates conductive heat transfer, which leads to moderate maximum surface temperatures. On the other hand, the high thermal diffusivity of aluminium results in an enlarged melt pool, see Figure 9, particularly in a larger melt pool depth, when compared with the melt pool dimensions predicted for stainless steel in Figure 4.

Furthermore, the maximum velocity magnitude trends determined for the aluminium substrate, see Figure 10, indicate a more distinct sensitivity to the laser fluence in comparison with the results for stainless steel presented in Figure 5. The maximum horizontal velocity magnitudes presented in Figure 10 for the aluminium melt, which notably exceed 10 m/s for elevated fluences, are at least twice as high as the ones in Figure 5 predicted for stainless steel. Contemplating on the Marangoni boundary condition in Equation (15), this difference is explained by the considerably smaller dynamic viscosity of liquid aluminium, which is only one fifth of the value for stainless steel according to Table 3. The higher velocity values for the liquid aluminium are even bounded by the lower magnitudes of the temperature coefficient of surface tension, cf. Table 3, and the horizontal temperature gradient, see Figure 11c. In conjunction with the markedly deeper melt pool, the increased horizontal melt velocities suggest that aluminium is a delicate substrate material, which requires a careful selection of the fluence for DLIP. Actually, the high melt velocities predicted in this work provide an explanation for the irregular surface morphology observed for aluminium substrates, see Figure 12c [96], unlike the nearly perfect repetitive topography of stainless steel, see Figure 12a [35], after DLIP at moderate fluences.

The surface microstructure and profile reprinted in Figure 12a,b were observed after a single pulse DLIP experiment with a periodicity of $\Lambda =5.05$ µm and a fluence of 0.6 J/cm^{2} on stainless steel. The width of the molten zone at the steel surface presented in Figure 12a is approximately 3.2 µm [35], which is in reasonable agreement with the maximum melt pool width of 2.9 µm, see Figure 4a, in the simulation. The surface profile in Figure 12b shows that melt is displaced from the interference maxima towards the interference minima [35]. The height difference between the original surface and the resulting depression is approximately 200 nm [35] and of the same order as the melt pool depth of 270 nm, see Figure 4b, computed in this work.

Concerning the high-purity aluminium substrate, the surface morphology and profile shown in Figure 12c,d, respectively, were obtained using a slightly smaller periodicity of $\Lambda =4.5$ µm and a fluence of 1.015 J/cm^{2}. For this reason, only the height of the resulting surface structures is compared with the melt pool depth predicted for aluminium. A structure height of approximately 0.8 µm can be inferred from the cross section in Figure 12d, which suggests that the volume of metal melt is considerably larger in the case of aluminium. The melt pool depth of 1.1 µm predicted in the present simulations for a fluence of 1.0 J/cm^{2}, see Figure 9b, is in agreement with the measured depth of the surface microstructure on aluminium after DLIP with a periodicity of $\Lambda =4.5$ µm, as reported in [96].

The present results are in reasonable agreement with previous work on DLIP of metals [12,96,97]. In detail, the maximum surface temperatures of aluminium during DLIP are in accordance with values computed using the FEM model, which were presented in [97]. However, lower maximum surface temperatures were predicted in [97] for stainless steel, probably owing to a considerably higher thermal diffusivity value employed. The magnitude of the horizontal temperature gradient at the surface of the stainless steel and aluminium substrates is in line with the thermal gradient determined in [97]. Further, the melt pool depth during DLIP of stainless steel presented in Figure 4b is of the same order of magnitude as the molten depth calculated in [12], in spite of the fact that the latter results were obtained for a smaller periodicity.

On the other hand, the molten depths determined from the thermal simulations for an aluminium substrate in [12] are notably lower than the present results in Figure 9b, although a comparable periodicity was considered. This deviation may be attributed to a higher thermal diffusivity of aluminium employed in [12], which entails a more homogeneous temperature field with less pronounced maxima and, therefore, a reduced absorption of the incident radiation. Nevertheless, the present depths of the aluminium melt pool are in the same size range as the structure depths measured after DLIP experiments on aluminium substrates using a similar periodicity, which were reported in [96].

Furthermore, the validity of the predicted velocity magnitudes is confirmed in the following. Considering the molten material as a shear layer, a characteristic horizontal velocity of the thermocapillary motion is obtained from the Marangoni boundary condition in Equation (15) as [98,99]
where ${\delta}_{\mathrm{m}}$ denotes a height related to the melt pool. If the height ${\delta}_{\mathrm{m}}$ is specified as half of the molten film thickness, ${u}_{\mathrm{tc}}$ in Equation (79) represents an average horizontal velocity of melt displacement [12,98]. For the present simulations, excessive horizontal velocity values would result if the melt pool depth was employed in Equation (79). Instead, it is suggested to substitute the thickness of the velocity boundary layer near the surface for ${\delta}_{\mathrm{m}}$. In addition, a time scale for thermocapillary convection is given by the ratio of half the interference pattern periodicity and the characteristic velocity magnitude [12,100]
where the horizontal temperature gradient is often assessed by relating a temperature difference to a width, e.g., $\partial T/\partial x\approx \Delta T/\Lambda /2$ [12], the resulting expression being proportional to ${\Lambda}^{2}/4$ [100].

$${u}_{\mathrm{tc}}=\frac{{\delta}_{\mathrm{m}}}{\eta}\frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}\gamma}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}T}\frac{\partial T}{\partial x},$$

$${t}_{\mathrm{tc}}=\frac{\Lambda /2}{\left|{u}_{\mathrm{tc}}\right|}=\frac{\eta \Lambda}{2{\delta}_{\mathrm{m}}\left|\frac{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}\gamma}{\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}T}\frac{\partial T}{\partial x}\right|},$$

Considering the simulations presented in detail in Section 5, the horizontal velocity magnitudes at the time $t=57$ ns are estimated by Equation (79) using the material properties from Table 3. For stainless steel subject to the laser fluence $2{\mathsf{\Phi}}_{0}=0.5$ J/cm^{2}, the average temperature gradient $\partial T/\partial x\approx -1240$ K/µm and ${\delta}_{\mathrm{m}}\approx 68$ nm result in ${u}_{\mathrm{tc}}\approx 5.13$ m/s. In case of the high-purity aluminium substrate and the energy density $2{\mathsf{\Phi}}_{0}=1.0$ J/cm^{2}, the average temperature gradient $\partial T/\partial x\approx -900$ K/µm and ${\delta}_{\mathrm{m}}\approx 115$ nm yield ${u}_{\mathrm{tc}}\approx 20.6$ m/s. The maximum surface velocity magnitudes computed in the ISPH simulations, as may be identified from Figure 6c and Figure 11c, are 32% and 15%, respectively, lower than these values. These deviations may be attributed to the transient character of the surface temperature distribution and the short duration of melt pool convection. Employing $\Lambda /2=2.5$ µm and the aforementioned velocities in Equation (80), the time scale for thermocapillary flow equals ${t}_{\mathrm{tc}}\approx 487$ ns for stainless steel and ${t}_{\mathrm{tc}}\approx 121$ ns for aluminium.

The present work demonstrates the use of SPH in the simulation of heat transfer and fluid flow during nanosecond pulsed DLIP of metallic substrates on the basis of a dimensionless formulation of the governing equations. For single-pulse treatment, the simulation results reveal a distinct behaviour of the dissimilar substrates stainless steel and high-purity aluminium. Particularly in the case of processing aluminium, the predictions confirm that thermocapillary convection is characterised by substantial velocity magnitudes beyond 10 m/s even at moderate laser fluences. Consequently, this outward flow from the centre of the melt pool surface, at the interference maximum, towards its edges is a conceivable structuring mechanism effective during DLIP of aluminium using a nanosecond pulse at low and moderate energy densities.

On the contrary, the higher absorptivity and lower thermal conductivity of stainless steel lead to high surface temperatures up to the vapourisation point near the interference maximum, even at moderate laser fluences. In addition, thermocapillary convection is less pronounced in liquid steel owing to its high dynamic viscosity. Therefore, it is more difficult to distinguish between the effects of the thermocapillary flow and the recoil pressure induced by vapourisation on the melt displacement during DLIP of stainless steel at moderate energy densities. The further investigation of the roles of thermocapillary convection and vapourisation-induced recoil pressure in the structuring of metal surfaces by DLIP ultimately necessitates the modelling of the melt pool surface deformation.

The numerical results presented in this work are compared with surface microstructures obtained by material characterisation after DLIP experiments on stainless steel and high-purity aluminium. In agreement with the observed surface morphologies, the simulations indicate a significantly deeper molten bath and a more effective melt displacement mechanism for the aluminium substrate. The predicted velocity magnitudes, which suggest a notable outward flow due to thermocapillary forces at the surface of the aluminium melt pool, are confirmed by a theoretical estimation.

In this work, an ISPH approach is deliberately employed to simulate the melt flow during DLIP of metals. This choice is motivated by virtue of the meaningful, physical pressure field computed using the ISPH technique. It is expected that an accurate pressure field is essential for the further study of melt displacement during DLIP, particularly for modelling the effect of the recoil pressure. The present model could serve as a basis for the prospective investigation. Although the SPH method is inherently capable of describing the evolution of a free surface in fluid flow, the application of a projection-based ISPH approach involves fundamental difficulties in this context. Provided that these issues can be eliminated, the method presented herein can be appropriately extended to simulate the melt displacement during DLIP of metals as well.

The research was conceptualised by A.F.L. and C.D., the methodology and software were developed by C.D., the results were discussed by A.F.L. and C.D., the original draft was written by C.D., the manuscript was reviewed and edited by A.F.L. and C.D. All authors have read and agreed to the published version of the manuscript.

This research received no external funding. The APC was borne in the frame of Open Access Funding by the Publication Fund of the TU Dresden and the Saxon State and University Library Dresden (SLUB).

The authors thank Achim Mahrle (TU Dresden, Institute of Manufacturing Technology) for the original suggestion of applying SPH to the modelling of laser material processing. Further, the first author acknowledges Theresa Jähnig (TU Dresden, Institute of Manufacturing Technology), Hartmut Krause, Subhashis Ray and Eric Werzner (TU Bergakademie Freiberg, Institute of Thermal Engineering) for the helpful discussions.

The authors declare no conflicts of interest.

The following abbreviations are used in this manuscript:

AISI | American Iron and Steel Institute |

BiCGStab | biconjugate gradient stabilised |

CG | conjugate gradient |

DLIP | direct laser interference patterning |

FEM | finite element method |

FWHM | full width at half maximum |

ISPH | incompressible smoothed particle hydrodynamics |

PPE | pressure Poisson equation |

SPH | smoothed particle hydrodynamics |

WCSPH | weakly compressible smoothed particle hydrodynamics |

- Berger, J.; Grosse Holthaus, M.; Pistillo, N.; Roch, T.; Rezwan, K.; Lasagni, A.F. Ultraviolet laser interference patterning of hydroxyapatite surfaces. Appl. Surf. Sci.
**2011**, 257, 3081–3087. [Google Scholar] [CrossRef] - Langheinrich, D.; Yslas, E.; Broglia, M.; Rivarola, V.; Acevedo, D.; Lasagni, A. Control of cell growth direction by direct fabrication of periodic micro- and submicrometer arrays on polymers. J. Polym. Sci., Part B Polym. Phys.
**2012**, 50, 415–422. [Google Scholar] [CrossRef] - Valle, J.; Burgui, S.; Langheinrich, D.; Gil, C.; Solano, C.; Toledo-Arana, A.; Helbig, R.; Lasagni, A.; Lasa, I. Evaluation of surface microtopography engineered by direct laser interference for bacterial anti-biofouling. Macromol. Biosci.
**2015**, 15, 1060–1069. [Google Scholar] [CrossRef] [PubMed] - Müller-Meskamp, L.; Kim, Y.H.; Roch, T.; Hofmann, S.; Scholz, R.; Eckhardt, S.; Leo, K.; Lasagni, A.F. Efficiency enhancement of organic solar cells by fabricating periodic surface textures using direct laser interference patterning. Adv. Mater.
**2012**, 24, 906–910. [Google Scholar] [CrossRef] [PubMed] - Roch, T.; Weihnacht, V.; Scheibe, H.J.; Roch, A.; Lasagni, A.F. Direct laser interference patterning of tetrahedral amorphous carbon films for tribological applications. Diamond Relat. Mater.
**2013**, 33, 20–26. [Google Scholar] [CrossRef] - Bieda, M.; Schmädicke, C.; Roch, T.; Lasagni, A. Ultra-low friction on 100Cr6-steel surfaces after direct laser interference patterning. Adv. Eng. Mater.
**2015**, 17, 102–108. [Google Scholar] [CrossRef] - Estevam-Alves, R.; Günther, D.; Dani, S.; Eckhardt, S.; Roch, T.; Mendonca, C.R.; Cestari, I.N.; Lasagni, A.F. UV direct laser interference patterning of polyurethane substrates as tool for tuning its surface wettability. Appl. Surf. Sci.
**2016**, 374, 222–228. [Google Scholar] [CrossRef] - Daniel, C.; Mücklich, F.; Liu, Z. Periodical micro-nano-structuring of metallic surfaces by interfering laser beams. Appl. Surf. Sci.
**2003**, 208–209, 317–321. [Google Scholar] [CrossRef] - Dowden, J.M. (Ed.) The Theory of Laser Materials Processing: Heat and Mass Transfer in Modern Technology, 1st ed.; Springer Series in Materials Science; Springer: Dordrecht, The Netherlands, 2009; Volume 119. [Google Scholar] [CrossRef]
- Daniel, C.; Lasagni, A.; Mücklich, F. Stress and texture evolution of Ni/Al multi-film by laser interference irradiation. Surf. Coat. Technol.
**2004**, 180–181, 478–482. [Google Scholar] [CrossRef] - Lasagni, A.; Mücklich, F. Study of the multilayer metallic films topography modified by laser interference irradiation. Appl. Surf. Sci.
**2005**, 240, 214–221. [Google Scholar] [CrossRef] - Lasagni, A.; D’Alessandria, M.; Giovanelli, R.; Mücklich, F. Advanced design of periodical architectures in bulk metals by means of Laser Interference Metallurgy. Appl. Surf. Sci.
**2007**, 254, 930–936. [Google Scholar] [CrossRef] - Lasagni, A.; Mücklich, F. FEM simulation of periodical local heating caused by laser interference metallurgy. J. Mater. Process. Technol.
**2009**, 209, 202–209. [Google Scholar] [CrossRef] - Gingold, R.A.; Monaghan, J.J. Smoothed particle hydrodynamics: Theory and application to non-spherical stars. Mon. Not. R. Astron. Soc.
**1977**, 181, 375–389. [Google Scholar] [CrossRef] - Lucy, L.B. A numerical approach to the testing of the fission hypothesis. Astron. J.
**1977**, 82, 1013–1024. [Google Scholar] [CrossRef] - Marongiu, J.C.; Leboeuf, F.; Caro, J.; Parkinson, E. Free surface flows simulations in Pelton turbines using an hybrid SPH-ALE method. J. Hydraul. Res.
**2010**, 48, 40–49. [Google Scholar] [CrossRef] - Pugliese Carratelli, E.; Viccione, G.; Bovolin, V. Free surface flow impact on a vertical wall: A numerical assessment. Theor. Comput. Fluid Dyn.
**2016**, 30, 403–414. [Google Scholar] [CrossRef] - Monaghan, J.J. Smoothed particle hydrodynamics and its diverse applications. Annu. Rev. Fluid Mech.
**2012**, 44, 323–346. [Google Scholar] [CrossRef] - Violeau, D. Fluid Mechanics and the SPH Method: Theory and Applications; Oxford University Press: Oxford, UK, 2012. [Google Scholar]
- Violeau, D.; Rogers, B.D. Smoothed particle hydrodynamics (SPH) for free-surface flows: Past, present and future. J. Hydraul. Res.
**2016**, 54, 1–26. [Google Scholar] [CrossRef] - Shadloo, M.S.; Oger, G.; Le Touzé, D. Smoothed particle hydrodynamics method for fluid flows, towards industrial applications: Motivations, current state, and challenges. Comput. Fluids
**2016**, 136, 11–34. [Google Scholar] [CrossRef] - Chen, J.; Beraun, J. Numerical study of ultrashort laser pulse interactions with metal films. Numer. Heat Transf. Part A
**2001**, 40, 1–20. [Google Scholar] - Gross, M. Smooth particle hydrodynamics (SPH) modelling of laser cutting. In International Congress on Applications of Lasers & Electro Optics; Laser Institute of America: Orlando, FL, USA, 2008; pp. 637–644. [Google Scholar]
- Muhammad, N.; Rogers, B.D.; Li, L. Understanding the behaviour of pulsed laser dry and wet micromachining processes by multi-phase smoothed particle hydrodynamics (SPH) modelling. J. Phys. D Appl. Phys.
**2013**, 46, 095101. [Google Scholar] [CrossRef] - Abidou, D.; Yusoff, N.; Nazri, N.; Awang, M.O.; Hassan, M.A.; Sarhan, A.A. Numerical simulation of metal removal in laser drilling using symmetric smoothed particle hydrodynamics. Precis. Eng.
**2017**, 49, 69–77. [Google Scholar] [CrossRef] - Tong, M.; Browne, D.J. Smoothed particle hydrodynamics modelling of the fluid flow and heat transfer in the weld pool during laser spot welding. IOP Conf. Ser. Mater. Sci. Eng.
**2012**, 27, 012080. [Google Scholar] [CrossRef] - Hu, H.; Eberhard, P. Thermomechanically coupled conduction mode laser welding simulations using smoothed particle hydrodynamics. Comput. Part. Mech.
**2017**, 4, 473–486. [Google Scholar] [CrossRef] - Russell, M.A.; Souto-Iglesias, A.; Zohdi, T.I. Numerical simulation of laser fusion additive manufacturing processes using the SPH method. Comput. Methods Appl. Mech. Eng.
**2018**, 341, 163–187. [Google Scholar] [CrossRef] - Yan, Y.; Li, L.; Sezer, K.; Wang, W.; Whitehead, D.; Ji, L.; Bao, Y.; Jiang, Y. CO
_{2}laser underwater machining of deep cavities in alumina. J. Eur. Ceram. Soc.**2011**, 31, 2793–2807. [Google Scholar] [CrossRef] - Hu, H.; Fetzer, F.; Berger, P.; Eberhard, P. Simulation of laser welding using advanced particle methods. GAMM Mitt. Ges. Angew. Math. Mech.
**2016**, 39, 149–169. [Google Scholar] [CrossRef] - Tanaka, M.; Cardoso, R.; Bahai, H. Modification of the LSMPS method for the conservation of the thermal energy in laser irradiation processes. Int. J. Numer. Methods Eng.
**2019**, 117, 161–187. [Google Scholar] [CrossRef] - Weirather, J.; Rozov, V.; Wille, M.; Schuler, P.; Seidel, C.; Adams, N.A.; Zaeh, M.F. A smoothed particle hydrodynamics model for laser beam melting of Ni-based alloy 718. Comput. Math. Appl.
**2019**, 78, 2377–2394. [Google Scholar] [CrossRef] - Liu, S.; Liu, J.; Chen, J.; Liu, X. Influence of surface tension on the molten pool morphology in laser melting. Int. J. Therm. Sci.
**2019**, 146, 106075. [Google Scholar] [CrossRef] - Fürstenau, J.P.; Wessels, H.; Weißenfels, C.; Wriggers, P. Generating virtual process maps of SLM using powder-scale SPH simulations. Comput. Part. Mech.
**2019**, 1–23. [Google Scholar] [CrossRef] - Demuth, C.; Bieda, M.; Lasagni, A.F.; Mahrle, A.; Wetzig, A.; Beyer, E. Thermal simulation of pulsed direct laser interference patterning of metallic substrates using the smoothed particle hydrodynamics approach. J. Mater. Process. Technol.
**2012**, 212, 689–699. [Google Scholar] [CrossRef] - Cao, Y.; Shin, Y.C. Multi-scale modeling of phase explosion in high fluence nanosecond laser ablation and clarification of ablation depth prediction criterion. Appl. Surf. Sci.
**2015**, 357, 74–85. [Google Scholar] [CrossRef] - Alshaer, A.W.; Rogers, B.D.; Li, L. Smoothed particle hydrodynamics (SPH) modelling of transient heat transfer in pulsed laser ablation of Al and associated free-surface problems. Comput. Mater. Sci.
**2017**, 127, 161–179. [Google Scholar] [CrossRef] - Komen, H.; Shigeta, M.; Tanaka, M. Numerical simulation of molten metal droplet transfer and weld pool convection during gas metal arc welding using incompressible smoothed particle hydrodynamics method. Int. J. Heat Mass Transf.
**2018**, 121, 978–985. [Google Scholar] [CrossRef] - Landau, L.D.; Lifshitz, E.M. Fluid Mechanics; Course of Theoretical Physics; Pergamon Press: London, UK, 1959; Volume 6. [Google Scholar]
- Nepomnyashchy, A.; Simanovskii, I.; Legros, J.C. Interfacial convection in Multilayer Systems; Springer: New York, NY, USA, 2006. [Google Scholar]
- Monaghan, J.J.; Lattanzio, J.C. A refined particle method for astrophysical problems. Astron. Astrophys.
**1985**, 149, 135–143. [Google Scholar] - Monaghan, J.J. An introduction to SPH. Comput. Phys. Commun.
**1988**, 48, 89–96. [Google Scholar] [CrossRef] - Price, D.J. Smoothed particle hydrodynamics and magnetohydrodynamics. J. Comput. Phys.
**2012**, 231, 759–794. [Google Scholar] [CrossRef] - Schoenberg, I.J. Contributions to the problem of approximation of equidistant data by analytic functions. Part A. On the problem of smoothing or graduation. A first class of analytic approximation formulae. Q. Appl. Math.
**1946**, 4, 45–99. [Google Scholar] [CrossRef] - Morris, J.P.; Fox, P.J.; Zhu, Y. Modeling low Reynolds number incompressible flows using SPH. J. Comput. Phys.
**1997**, 136, 214–226. [Google Scholar] [CrossRef] - Speith, R. Untersuchung von Smoothed Particle Hydrodynamics Anhand Astrophysikalischer Beispiele. Ph.D. Thesis, Eberhard Karls Universität Tübingen, Tübingen, Germany, 1998. (In German). [Google Scholar]
- Leroy, A. A New Incompressible SPH Model: Towards Industrial Applications (Un Nouveau modèle SPH Incompressible: Vers L’application à des cas Industriels). Ph.D. Thesis, Université Paris-Est, Paris, France, 2014. [Google Scholar]
- Monaghan, J.J. Smoothed particle hydrodynamics. Annu. Rev. Astron. Astrophys.
**1992**, 30, 543–574. [Google Scholar] [CrossRef] - Monaghan, J.J. Smoothed particle hydrodynamics. Rep. Prog. Phys.
**2005**, 68, 1703–1759. [Google Scholar] [CrossRef] - Cummins, S.J.; Rudman, M. An SPH projection method. J. Comput. Phys.
**1999**, 152, 584–607. [Google Scholar] [CrossRef] - Brookshaw, L. A method of calculating radiative heat diffusion in particle simulations. Proc. Astron. Soc. Aust.
**1985**, 6, 207–210. [Google Scholar] [CrossRef] - Cleary, P.W.; Monaghan, J.J. Conduction modelling using smoothed particle hydrodynamics. J. Comput. Phys.
**1999**, 148, 227–264. [Google Scholar] [CrossRef] - Randles, P.; Libersky, L.D. Normalized SPH with stress points. Int. J. Numer. Methods Eng.
**2000**, 48, 1445–1462. [Google Scholar] [CrossRef] - Randles, P.; Libersky, L.D. Smoothed particle hydrodynamics: Some recent improvements and applications. Comput. Methods Appl. Mech. Eng.
**1996**, 139, 375–408. [Google Scholar] [CrossRef] - Libersky, L.D.; Petschek, A.G. Smooth particle hydrodynamics with strength of materials. In Advances in the Free-Lagrange Method Including Contributions on Adaptive Gridding and the Smooth Particle Hydrodynamics Method; Springer: Berlin/Heidelberg, Germany, 1991; pp. 248–257. [Google Scholar]
- Monaghan, J.J. Simulating free surface flows with SPH. J. Comput. Phys.
**1994**, 110, 399–406. [Google Scholar] [CrossRef] - Kulasegaram, S.; Bonet, J.; Lewis, R.; Profit, M. A variational formulation based contact algorithm for rigid boundaries in two-dimensional SPH applications. Comput. Mech.
**2004**, 33, 316–325. [Google Scholar] [CrossRef] - Shao, S.; Lo, E.Y. Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface. Adv. Water Resour.
**2003**, 26, 787–800. [Google Scholar] [CrossRef] - Lee, E.S.; Moulinec, C.; Xu, R.; Violeau, D.; Laurence, D.; Stansby, P. Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method. J. Comput. Phys.
**2008**, 227, 8417–8436. [Google Scholar] [CrossRef] - Colagrossi, A.; Bouscasse, B.; Antuono, M.; Marrone, S. Particle packing algorithm for SPH schemes. Comput. Phys. Commun.
**2012**, 183, 1641–1653. [Google Scholar] [CrossRef] - Molteni, D.; Colagrossi, A. A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH. Comput. Phys. Commun.
**2009**, 180, 861–872. [Google Scholar] [CrossRef] - Antuono, M.; Colagrossi, A.; Marrone, S.; Molteni, D. Free-surface flows solved by means of SPH schemes with numerical diffusive terms. Comput. Phys. Commun.
**2010**, 181, 532–549. [Google Scholar] [CrossRef] - Sun, P.N.; Colagrossi, A.; Marrone, S.; Antuono, M.; Zhang, A.M. A consistent approach to particle shifting in the δ-Plus-SPH model. Comput. Methods Appl. Mech. Eng.
**2019**, 348, 912–934. [Google Scholar] [CrossRef] - Chorin, A.J. Numerical solution of the Navier-Stokes equations. Math. Comput.
**1968**, 22, 745–762. [Google Scholar] [CrossRef] - Temam, R. Une méthode d’approximation de la solution des équations de Navier-Stokes. Bull. Soc. Math. Fr.
**1968**, 96, 115–152. (In French) [Google Scholar] [CrossRef] - Chorin, A.J.; Marsden, J.E. A Mathematical Introduction to Fluid Mechanics, 3rd ed.; Texts in Applied Mathematics; Springer: New York, NY, USA, 1997; Volume 4. [Google Scholar]
- Hu, X.; Adams, N.A. An incompressible multi-phase SPH method. J. Comput. Phys.
**2007**, 227, 264–278. [Google Scholar] [CrossRef] - Ellero, M.; Serrano, M.; Español, P. Incompressible smoothed particle hydrodynamics. J. Comput. Phys.
**2007**, 226, 1731–1752. [Google Scholar] [CrossRef] - Xu, R.; Stansby, P.; Laurence, D. Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach. J. Comput. Phys.
**2009**, 228, 6703–6725. [Google Scholar] [CrossRef] - Allen, M.; Tildesley, D. Computer Simulation of Liquids; Oxford University Press: Oxford, UK, 1987. [Google Scholar]
- Hockney, R.; Eastwood, J. Computer Simulation Using Particles; McGraw–Hill: New York, NY, USA, 1981. [Google Scholar]
- Monaghan, J.J.; Gingold, R.A. Shock simulation by the particle method SPH. J. Comput. Phys.
**1983**, 52, 374–389. [Google Scholar] [CrossRef] - Viccione, G.; Bovolin, V.; Pugliese Carratelli, E. Defining and optimizing algorithms for neighbouring particle identification in SPH fluid simulations. Int. J. Numer. Methods Fluids
**2008**, 58, 625–638. [Google Scholar] [CrossRef] - Domínguez, J.M.; Crespo, A.J.C.; Gómez-Gesteira, M.; Marongiu, J.C. Neighbour lists in smoothed particle hydrodynamics. Int. J. Numer. Methods Fluids
**2011**, 67, 2026–2042. [Google Scholar] [CrossRef] - Morris, J.P. Analysis of Smoothed Particle Hydrodynamics with Applications. Ph.D. Thesis, Monash University, Melbourne, Australia, 1996. [Google Scholar]
- Hackbusch, W. Elliptic differential equations: Theory and numerical treatment; Springer Series in Computational Mathematics; Springer: Berlin/Heidelberg, Germany, 1992; Volume 18. [Google Scholar]
- Hestenes, M.R.; Stiefel, E. Methods of conjugate gradients for solving linear systems. J. Res. Natl. Bur. Stand.
**1952**, 49, 409–436. [Google Scholar] [CrossRef] - Saad, Y. Iterative methods for sparse linear systems, 2nd ed.; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2003. [Google Scholar]
- Hosseini, S.M.; Feng, J.J. Pressure boundary conditions for computing incompressible flows with SPH. J. Comput. Phys.
**2011**, 230, 7473–7487. [Google Scholar] [CrossRef] - van der Vorst, H.A. Bi-CGStab: A fast and smoothly convergent variant of Bi-CG for the solution of non-symmetric linear systems. SIAM J. Sci. Statist. Comput.
**1992**, 13, 631–644. [Google Scholar] [CrossRef] - Graves-Morris, P.R. The breakdowns of BiCGStab. Numer. Algorithms
**2002**, 29, 97–105. [Google Scholar] [CrossRef] - Sahoo, P.; Debroy, T.; McNallan, M.J. Surface tension of binary metal—Surface active solute systems under conditions relevant to welding metallurgy. Metall. Trans. B
**1988**, 19, 483–491. [Google Scholar] [CrossRef] - Mills, K.C. Recommended values of thermophysical properties for selected commercial alloys; Woodhead Publishing: Cambridge, UK, 2002. [Google Scholar]
- Steen, W.M.; Mazumder, J. Laser Material Processing, 4th ed.; Springer: London, UK, 2010. [Google Scholar]
- Lide, D.R. (Ed.) Handbook of Chemistry and Physics, 85th ed.; CRC Press: Boca Raton, FL, USA, 2004. [Google Scholar]
- Grigoriev, I.S.; Meilikhov, E.Z. (Eds.) Handbook of Physical Quantities; CRC Press: Boca Raton, FL, USA, 1997; Chapter 12; pp. 361–389. [Google Scholar]
- Gąsior, W.; Moser, Z.; Pstruś, J. Densities of solid aluminum-lithium (Al-Li) alloys. J. Phase Equilib.
**1998**, 19, 234–238. [Google Scholar] [CrossRef] - Assael, M.J.; Kakosimos, K.; Banish, R.M.; Brillo, J.; Egry, I.; Brooks, R.; Quested, P.N.; Mills, K.C.; Nagashima, A.; Sato, Y.; et al. Reference data for the density and viscosity of liquid aluminum and liquid iron. J. Phys. Chem. Ref. Data
**2006**, 35, 285–300. [Google Scholar] [CrossRef] - Valencia, J.J.; Quested, P.N. ASM Handbook, Vol. 15: Casting; Thermophysical Properties; ASM International: Materials Park, OH, USA, 2008; pp. 468–481. [Google Scholar]
- Touloukian, Y.S.; Powell, R.W.; Ho, C.Y.; Klemens, P.G. Thermal Conductivity—Metallic Elements and Alloys; Thermophysical Properties of Matter—The TPRC Data Series; IFI/Plenum: New York, NY, USA, 1970; Volume 1. [Google Scholar]
- Rothwell, E. A precise determination of the viscosity of liquid tin, lead, bismuth, and aluminium by an absolute method. J. Inst. Met.
**1962**, 90, 389–394. [Google Scholar] - Kaptay, G. A unified model for the cohesive enthalpy, critical temperature, surface tension and volumetric thermal expansion coefficient of liquid metals of bcc, fcc and hcp crystals. Mater. Sci. Eng. A
**2008**, 495, 19–26. [Google Scholar] [CrossRef] - Nasch, P.M.; Steinemann, S.G. Density and thermal expansion of molten manganese, iron, nickel, copper, aluminum and tin by means of the gamma-ray attenuation technique. Phys. Chem. Liq.
**1995**, 29, 43–58. [Google Scholar] [CrossRef] - Sarou-Kanian, V.; Millot, F.; Rifflet, J.C. Surface tension and density of oxygen-free liquid aluminum at high temperature. Int. J. Thermophys.
**2003**, 24, 277–286. [Google Scholar] [CrossRef] - Modest, M.F. Reflectivity and absorptivity of opaque surfaces. In Handbook of Laser Materials Processing; Ready, J.F., Farson, D.F., Eds.; Magnolia Publishing, Laser Institute of America: Orlando, FL, USA, 2001; pp. 175–182. [Google Scholar]
- D’Alessandria, M.; Lasagni, A.; Mücklich, F. Direct micropatterning of aluminum substrates via laser interference metallurgy. Appl. Surf. Sci.
**2008**, 255, 3210–3216. [Google Scholar] [CrossRef] - Bieda, M.; Beyer, E.; Lasagni, A.F. Direct fabrication of hierarchical microstructures on metals by means of direct laser interference patterning. J. Eng. Mater. Technol.
**2010**, 132, 031015. [Google Scholar] [CrossRef] - Cazabat, A.M.; Heslot, F.; Troian, S.M.; Carles, P. Fingering instability of spreading thin films driven by temperature gradients. Nature
**1990**, 346, 824–846. [Google Scholar] [CrossRef] - Bäuerle, D. Laser Processing and Chemistry, 4th ed.; Springer: Berlin/Heidelberg, Germany, 2011. [Google Scholar]
- Favazza, C.; Trice, J.; Kalyanaraman, R.; Sureshkumar, R. Self-organized metal nanostructures through laser-interference driven thermocapillary convection. Appl. Phys. Lett.
**2007**, 91, 043105. [Google Scholar] [CrossRef]

Particle | Role of Particle | |
---|---|---|

Fluid Particle | Edge Particle | |

i | $a}_{ii}=2\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^{{N}_{i}}\frac{{m}_{j}^{*}}{{\rho}_{j}^{*}}\frac{1}{{\overline{r}}_{ij}^{*}}\frac{\partial W\left({\overline{r}}_{ij}^{*},{l}_{\mathrm{sm}}^{*}\right)}{\partial {\overline{r}}_{ij}^{*}$ | $a}_{ii}=2\sum _{\begin{array}{c}j=1,j\ne i\\ j\notin \left\{{i}_{\mathrm{D}}\right\}\end{array}}^{{N}_{i}}\frac{{m}_{j}^{*}}{{\rho}_{j}^{*}}\frac{1}{{\overline{r}}_{ij}^{*}}\frac{\partial W\left({\overline{r}}_{ij}^{*},{l}_{\mathrm{sm}}^{*}\right)}{\partial {\overline{r}}_{ij}^{*}$ |

j | $a}_{ij}=-2\frac{{m}_{j}^{*}}{{\rho}_{j}^{*}}\frac{1}{{\overline{r}}_{ij}^{*}}\frac{\partial W\left({\overline{r}}_{ij}^{*},{l}_{\mathrm{sm}}^{*}\right)}{\partial {\overline{r}}_{ij}^{*}$ | $a}_{ij}=-2\frac{{m}_{j}^{*}}{{\rho}_{j}^{*}}\frac{1}{{\overline{r}}_{ij}^{*}}\frac{\partial W\left({\overline{r}}_{ij}^{*},{l}_{\mathrm{sm}}^{*}\right)}{\partial {\overline{r}}_{ij}^{*}}-2\sum _{{j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}\in \left\{{j}_{\mathrm{D}}\right\}}\frac{{m}_{{j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}}^{*}}{{\rho}_{{j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}}^{*}}\frac{1}{{\overline{r}}_{i{j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}}^{*}}\frac{\partial W\left({\overline{r}}_{i{j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}}^{*},{l}_{\mathrm{sm}}^{*}\right)}{\partial {\overline{r}}_{i{j}^{\prime}\phantom{\rule{-0.66666pt}{0ex}}}^{*}$ |

Process Parameter | Symbol | Value |
---|---|---|

wavelength | $\lambda $ | 355 nm |

full angle between beams | $\theta $ | 0.071 rad |

periodicity of interference pattern | $\Lambda $ | 5 µm |

energy density (fluence) per beam | ${\Phi}_{0}$ | 0.3 J/cm^{2} |

pulse duration (FWHM) | ${\tau}_{\mathrm{p}}$ | 10 ns |

pulse time | ${t}_{\mathrm{p}}$ | 50 ns |

simulation duration | ${t}_{\mathrm{end}}$ | 200 ns |

initial substrate temperature | ${T}_{0}$ | 298.15 K |

gravitational acceleration | g | 9.81 m/s^{2} |

Fourier number | $\mathit{Fo}$ | 5.0 |

Material Property | Symbol | AISI 304 Stainless Steel | High-Purity Aluminium | Unit | Refs. |
---|---|---|---|---|---|

solidus temperature | ${T}_{\mathrm{s}}$ | 1673 | ${T}_{\mathrm{m}}$: 933.35 | K | [83] |

liquidus temperature | ${T}_{\mathrm{l}}$ | 1727 | K | [83] | |

vapourisation temperature | ${T}_{\mathrm{v}}$ | 3273 | 2792 | K | [84,85] |

enthalpy of fusion | ${L}_{\mathrm{f}}$ | 251 | 397 | kJ/kg | [83,84] |

enthalpy of vapourisation | ${L}_{\mathrm{v}}$ | 6500 | 10860 | kJ/kg | [84,86] |

density | ${\rho}_{\mathrm{ref}}$ | 7262 | 2228 | kg/m^{3} | [83,87,88] |

specific heat | ${c}_{\mathrm{p},\mathrm{ref}}$ | 704 | 1077 | J/(kgK) | [83,89] |

thermal conductivity | ${\kappa}_{\mathrm{ref}}$ | 26.8 | 139.5 | W/(mK) | [83,90] |

thermal diffusivity | ${a}_{\mathrm{ref}}$ | $5.24\xb7{10}^{-6}$ | $5.75\xb7{10}^{-5}$ | m^{2}/s | |

dynamic viscosity (at ${T}_{\mathrm{l}}$ or ${T}_{\mathrm{m}}$) | $\eta $ | $7.0\xb7{10}^{-3}$ | $1.38\xb7{10}^{-3}$ | Pas | [83,91] |

kinematic viscosity (at ${T}_{\mathrm{l}}$ or ${T}_{\mathrm{m}}$) | $\nu $ | $1.02\xb7{10}^{-6}$ | $5.80\xb7{10}^{-7}$ | m^{2}/s | |

volumetric thermal expansion coefficient | $\beta $ | $8.5\xb7{10}^{-5}$ | $9.8\xb7{10}^{-5}$ | 1/K | [92,93] |

temperature coefficient of surface tension | $\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}\gamma /\phantom{\rule{-0.166667em}{0ex}}\mathrm{d}T$ | $-4.3\xb7{10}^{-4}$ | $-2.75\xb7{10}^{-4}$ | N/(mK) | [82,94] |

absorption coefficient (at 355 nm) | $\alpha $ | $9.56\xb7{10}^{7}$ | $1.52\xb7{10}^{8}$ | 1/m | [85] |

reflectivity (at 355 nm) | R | $0.556$ | $0.925$ (at ${T}_{0}$) | – | [85,95] |

Quantity | Symbol | AISI 304 Stainless Steel | High-Purity Aluminium |
---|---|---|---|

thermal diffusion length | L | 457.8 nm | 1.517 µm |

Laser number | $\mathit{La}$ | 37.7151 | 152.657 |

solid-liquid phase change number | ${\mathit{Ph}}_{\mathrm{s}/\mathrm{l}}$ | 0.119849 | 0.147849 |

liquid-vapour phase change number | ${\mathit{Ph}}_{\mathrm{l}/\mathrm{v}}$ | 3.10367 | 4.04443 |

Prandtl number | $\mathit{Pr}$ | 0.1947 | 0.0101 |

Rayleigh number | $\mathit{Ra}$ | $2.31448\xb7{10}^{-8}$ | $1.86969\xb7{10}^{-7}$ |

Marangoni number | $\mathit{Ma}$ | 8.29744 | 9.76488 |

Substrate | Quantity | Value | Unit | ||||||
---|---|---|---|---|---|---|---|---|---|

AISI 304 | $2{\Phi}_{0}$ | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 | J/cm^{2} |

$\mathit{La}$ | 25.1434 | 31.4292 | 37.7151 | 44.0009 | 50.2868 | 56.5726 | 62.8585 | – | |

fine zone | 270 | 270 | 320 | 320 | 320 | 370 | 370 | nm | |

Al 99.99% | $2{\Phi}_{0}$ | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 | 1.1 | 1.2 | J/cm^{2} |

$\mathit{La}$ | 152.657 | 178.100 | 203.543 | 228.986 | 254.429 | 279.871 | 305.314 | – | |

fine zone | 0.77 | 0.77 | 1.02 | 1.02 | 1.27 | 1.27 | 1.52 | µm |

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).