**2.2. Laser applications and their simulation in the literature**

#### *2.2.1. Laser welding*

process, numerical simulation can give an insight of what is happening inside the work‐

Multiphysics packages such as ANSYS, ABAQUS, ADINA, or COMSOL are effective for the aforementioned issue and for fundamental investigation. In the industry, another semiempir‐ ical modeling method is useful when direct answers are expected regarding parameters optimization, predicting the system behavior or analyzing the effect of a modification. The Design Of Experiment (DOE) methodology shows its efficiency especially when the experi‐ mental investigation suffers from constraints such as the availability of the device, the

Although nonexhaustive, a list of laser applications and how they were modeled will be presented in the following parts. Multiphysics models will be presented first, followed by the DOE approach. A laser application will be discussed for both types of numerical simulations.

In general, for a solid, the laser ablation process strongly depends on the absorption of photons by the material. Hence, it relies on the laser wavelength, its pulse duration, its fluence (energy density per surface area), and also on the material properties. The photon energy is absorbed and converted into kinetic energy by vibration of the electron cloud. This energy is then transmitted to the volume by phonon-electron coupling [7, 8]. Above picosecond pulse duration, a laser beam can be treated as a heat source. Hence, whether it is drilling, cutting, welding…, a laser process can be seen as a thermodynamic problem and the classical heat transfer equation (Eq. (4)) can be used to determine the spatial and temporal distribution of the temperature in the workpiece. Then for advanced studies, the Navier-Stokes equations can be implemented when phase transition is involved (Eq. (2)). The final thermo-mechanical state of the substrate or the plasma generation during laser ablation and its effect on the efficiency

To be precise, there should be one equation describing the temperature evolution for the electron network and another one for the ions lattice (Section 2.3.1). Indeed, the characteristic time of energy transfer between electrons and ions is between 0.1 a few ps. For pulse duration longer than 500 fs, electrons and ions are both heated during the laser pulse. Hence, one equation is sufficient. Below 500 fs, the lattice is poorly affected because it is decoupled with the electrons and two equations are required [9]. However, as seen in the literature, one temperature model could be used to describe glass welding [10] or the formation of wave‐ guides with femtosecond laser when high frequency is used (hundreds of kHz to MHz range) [11]. After a short review of the use of multiphysics simulation, the heat transfer model will be detailed and will lead to the investigation of picosecond laser induced periodic surface

piece, which is delicate by experimental means.

162 Modeling and Simulation in Engineering Sciences

experimenter, the material, the risks, the costs, or the environment.

**2. Multiphysics simulation of laser processing**

of the laser process can be investigated as well.

structuring (LIPSS) on copper film.

**2.1. Introduction**

One of the main processes to be modeled is probably laser welding since it is widely used in various industries. Laser is contactless, repeatable, automated, generates low distortion due to heat, and generates high throughput. The main concern for the manufacturer is the quality of the weld, which is characterized by the geometry of the bead (the bead width and the depth of penetration of the weld) and also the tensile strength of the joint. The bead geometry is directly linked to the laser parameters which will determine how and what quantity of material is melted. Moreover, the heating and cooling cycle can be analyzed and this is relevant for the study of residual stresses and more importantly for the mechanical strength of the joint.

This kind of analysis is performed by Balasubramanian et al. [12], for instance, during the study of laser welding of AISI304 stainless steel sheet with the finite element method (FEM) package SYSWELD. A 3D conical Gaussian beam was modeled, thermal dependency of material properties was taken into account, and the laser was a solid-state Nd:YAG operating in CW (continuous wave) mode. The important part was about the Gaussian beam, which had to be modeled according a specific way to take into account the high depth/width aspect-ratio of a welding process. Former models would use a simple coefficient 0<*A*<1 for the surface absorp‐ tivity and it would underestimate the penetration depth of the weld because it is a volume process and not a surface process. Also, the thermal property dependency with temperature is important for the thermal diffusion and therefore the final analysis of the HAZ and the weld depth. Satisfactory results were obtained since the standard deviation between the experi‐ mental weld pool geometry and the simulation one was 4.54%. The laser parameters were a laser power of 1250 W, a scanning speed of 750 mm/min, and a beam angle of 90° (**Figure 1**).

**Figure 1.** Comparison of fusion zone. (a) Microstructure image, (b) FEM analysis [12].

#### *2.2.2. Laser ablating process*

Andreas Otto et al. went a step further with the integration of fluid dynamics [1]. They wanted to discuss the dynamics of the process such as welding or evaporative cutting (keyhole/melt pool oscillations and vapor flows) with the help of simulation. They combined the use of OpenFOAM (Open Field Operation and Manipulation) and COMSOL simulation to do so. OpenFOAM was used to solve the fluid dynamics part by the means of Volume Of Fluid (VOF) approach. Automatic remeshing is also implemented to cope up with material deformation during melting and evaporation. The VOF method is used to track the deforming surfaces and open boundaries. A fraction function is defined and is set at zero when the computational cell is only solid and increases (up to 1) as the cell is filled with liquid. In this chapter, the VOF method was adjusted so that vapor flow could be evaluated as well. Indeed, a compression algorithm for cells marked as vapor is activated by default in the software and causes discon‐ tinuities in the fraction function and thus in the resulting tracing of the interfaces. The position of the boundary relative to the mesh can be tracked by the Level Set method. This method also uses a tracking variable, but on the entire domain this time. However, VOF is preferred as the mass of the traced cell is conserved. The electromagnetic wave propagation (laser beam propagating inside the material and potential vapor formed during laser ablation) and thermomechanical part were simulated with COMSOL (finite element method, FEM). It is also used to couple the equations.

The fluid dynamics is modeled according to the Navier-Stokes equation (Eq. (2)). The fluid is assumed to be incompressible as its speed does not exceed half of the sound velocity. However, vapor is compressible and a mixed model is planned for future investigation. Mass conserva‐ tion is assumed (Eq. (1)).

$$\frac{\partial \rho}{\partial t} + \nabla \rho \vec{u} = 0 \text{ (mass conservation equation)}\tag{1}$$

$$
\rho \frac{\partial \vec{u}}{\partial t} + \rho \vec{u} \nabla \vec{u} = \cdot \rho \nabla \mathbf{P} + \eta \Delta \vec{u} \text{ (Navier-Stokes)}\tag{2}
$$

$$P\_0 + P\_\mathbf{v} + P\_\mathbf{L} = 2\mathbf{K}\sigma\tag{3}$$

where ρ is the material density, *t* is the time variable, *u* is the velocity field, *P*0 is the atmospheric pressure, *P*v is the vapor pressure, *P*L is the pressure in the liquid that interacts with the surface tension σ caused by the curvature *K* of the metal surface. The surface tension depends on the fluid temperature and is also important for tracking the interface deformation.

The temperature appears in the vapor pressure equation and thus in the Navier-Stokes equation (Eq. (2)). It needs to be coupled with the heat transfer equation (Eq. (4)) that will describe the temperature variation during laser irradiation. The temperature variation will have an effect on the fluid flow as understood from Navier-heat transfer equation coupling. As mentioned before, in this case, the laser beam can be modeled as a heat source *Q*:

*2.2.2. Laser ablating process*

164 Modeling and Simulation in Engineering Sciences

to couple the equations.

tion is assumed (Eq. (1)).

*t* r

¶

r

¶

r

 r

¶ +Ñ =

Andreas Otto et al. went a step further with the integration of fluid dynamics [1]. They wanted to discuss the dynamics of the process such as welding or evaporative cutting (keyhole/melt pool oscillations and vapor flows) with the help of simulation. They combined the use of OpenFOAM (Open Field Operation and Manipulation) and COMSOL simulation to do so. OpenFOAM was used to solve the fluid dynamics part by the means of Volume Of Fluid (VOF) approach. Automatic remeshing is also implemented to cope up with material deformation during melting and evaporation. The VOF method is used to track the deforming surfaces and open boundaries. A fraction function is defined and is set at zero when the computational cell is only solid and increases (up to 1) as the cell is filled with liquid. In this chapter, the VOF method was adjusted so that vapor flow could be evaluated as well. Indeed, a compression algorithm for cells marked as vapor is activated by default in the software and causes discon‐ tinuities in the fraction function and thus in the resulting tracing of the interfaces. The position of the boundary relative to the mesh can be tracked by the Level Set method. This method also uses a tracking variable, but on the entire domain this time. However, VOF is preferred as the mass of the traced cell is conserved. The electromagnetic wave propagation (laser beam propagating inside the material and potential vapor formed during laser ablation) and thermomechanical part were simulated with COMSOL (finite element method, FEM). It is also used

The fluid dynamics is modeled according to the Navier-Stokes equation (Eq. (2)). The fluid is assumed to be incompressible as its speed does not exceed half of the sound velocity. However, vapor is compressible and a mixed model is planned for future investigation. Mass conserva‐

*u* 0 (mass conservation equation)

s

where ρ is the material density, *t* is the time variable, *u* is the velocity field, *P*0 is the atmospheric pressure, *P*v is the vapor pressure, *P*L is the pressure in the liquid that interacts with the surface tension σ caused by the curvature *K* of the metal surface. The surface tension depends on the

The temperature appears in the vapor pressure equation and thus in the Navier-Stokes equation (Eq. (2)). It needs to be coupled with the heat transfer equation (Eq. (4)) that will describe the temperature variation during laser irradiation. The temperature variation will


0vL *PPP* + + =K2

fluid temperature and is also important for tracking the interface deformation.

 rh

¶ + Ñ = Ñ+D

<sup>r</sup> (1)

(3)

<sup>r</sup> rr r (2)

$$\frac{\partial \rho H}{\partial t} - \nabla k \nabla T + \nabla \rho \vec{u} H = \underline{Q} \tag{4}$$

$$H = L + C\_p \times T \tag{5}$$

where *T* is the temperature, *C*p is the heat capacity, *k* is the heat conductivity, *H* is the total enthalpy, and *L* is the latent heat of phase change. *L* is adjusted step by step by comparing the new result with the previous one until the difference is lower than a defined threshold (it takes approximately 5–10 iteration steps with *L* being neglected at first). Eq. (5) is thus used to evaluate the mass flow which is implemented as an input in the Navier-Stokes equation to simulate the fluid flow.

Important observations were made. For a deep welding process for instance, it is known that the melt pool will oscillate and cause keyhole width variation and porosity might appear depending on laser parameters and render the joint more fragile. The simulation reveals that waves of molten material run down the front of the keyhole. They cause a periodical modifi‐ cation of the keyhole diameter (keyhole oscillations). The liquid is accelerated around the keyhole. At about two thirds of the length of the melt pool, the liquid hits the backflow from the back of the melt pool and turbulences are observed in the lower rear part. On the other hand, the flow pattern is more laminar in the upper part of the melt pool. It causes the melt flow to change direction and reclose the keyhole. At higher scanning speeds, the simulation reveals the formation of pores.

For ablative processes (**Figure 2**) such as cutting, drilling, or structuring, it was found that droplets would be expulsed from the groove in a periodical way while using laser intensities above 108 W/cm2 .

The same authors studied the effect of the laser wavelength on the dynamics of a welding process [13]. Indeed, the temperature variation of the material depends on the absorption of the laser beam. The propagation of the beam also depends on the topology changes during melting and evaporation. Changes of propagation and absorption will impact the attenuation or amplification of the keyhole instabilities discussed previously. The pore formation (and thus the weld quality) depends on these oscillations and on the recoil pressure exerted on the weld pool when the material is evaporated. In order to study the beam absorption effect on the welding quality, they used a CO2 laser operating at a wavelength of 10.6 μm and a Nd:YAG laser operating at 1.064 μm. The feed rate was fixed at 0.1 mm/s, the laser power at 4.2 kW, and the beam radius at 200 μm. Compared to the previous model, the heat source term was modified so that Fresnel equations were taken into account. They took into account the absorption dependency on wavelength, angle of incidence, and polarization of the beam. It was found that oscillations are strongly influenced by fluctuations at the keyhole front for a wavelength of 1.064 μm. At 10.6 μm, most of the laser energy is absorbed at the keyhole front, provided that the surface of the molten steel is parallel to the incident laser beam. The absorption is reduced when parts of the keyhole front are nearly normal to the incident laser beam.

**Figure 2.** Simulation of laser beam drilling with short pulsed lasers [1].

The same kind of model coupling the heat transfer and the Navier-Stokes equations is used by other authors [14, 15]. In Ref. [14], the convection induced by the surface state and volume force is represented using the Boussinesq approximation (which states that the fluid is partially incompressible). The Marangoni effect (surface tension gradient depends on temperature gradient) was also taken into account. In Ref. [15], the Marangoni effect is neglected and the melt pool is considered as an incompressible fluid. The Level Set method is chosen as it is said to be efficient in treating the complex changes at the interfaces. Thus, it can model more easily the pores formation in the weld bead. Another approach was used to predict the ablated crater geometry in a PMMA substrate. Radice et al. [16] used a phase change tracking equation *h*(*T*) (inspired from the VOF method) in COMSOL Multiphysics to monitor the material modifica‐ tion over the phase change temperature range (Eq. (6)).

$$h(T) = f \text{lc2} \text{hs}(T - T1, T2 - T1) \tag{6}$$

where *T*1 is the temperature where the phase change begins, *T*2 is the temperature where it ends, and flc2hs is a smoothed Heaviside function used in COMSOL. *T*1, *T*2, and the latent heat of phase change were determined by a differential scanning calorimeter (DSC). *h*(*T*) is null when the material is in its solid state and it equals 1 when the material vaporizes. Coupled with the conventional heat transfer equation (latent heat of phase change is taken into account) (Eq. (4)), the authors were able to determine the hole geometry which revealed to be in good agreement with experimental results.

## *2.2.3. Example of nonablative process*

provided that the surface of the molten steel is parallel to the incident laser beam. The absorption is reduced when parts of the keyhole front are nearly normal to the incident laser

The same kind of model coupling the heat transfer and the Navier-Stokes equations is used by other authors [14, 15]. In Ref. [14], the convection induced by the surface state and volume force is represented using the Boussinesq approximation (which states that the fluid is partially incompressible). The Marangoni effect (surface tension gradient depends on temperature gradient) was also taken into account. In Ref. [15], the Marangoni effect is neglected and the melt pool is considered as an incompressible fluid. The Level Set method is chosen as it is said to be efficient in treating the complex changes at the interfaces. Thus, it can model more easily the pores formation in the weld bead. Another approach was used to predict the ablated crater geometry in a PMMA substrate. Radice et al. [16] used a phase change tracking equation *h*(*T*) (inspired from the VOF method) in COMSOL Multiphysics to monitor the material modifica‐

where *T*1 is the temperature where the phase change begins, *T*2 is the temperature where it ends, and flc2hs is a smoothed Heaviside function used in COMSOL. *T*1, *T*2, and the latent heat of phase change were determined by a differential scanning calorimeter (DSC). *h*(*T*) is null when the material is in its solid state and it equals 1 when the material vaporizes. Coupled with the conventional heat transfer equation (latent heat of phase change is taken into account) (Eq. (4)), the authors were able to determine the hole geometry which revealed to be in good

*h T flc2hs T T T T* ( ) ( 1, 2 1) = -- (6)

**Figure 2.** Simulation of laser beam drilling with short pulsed lasers [1].

tion over the phase change temperature range (Eq. (6)).

agreement with experimental results.

beam.

166 Modeling and Simulation in Engineering Sciences

Other applications where thermo-mechanical studies are important were also investigated. Glass cutting, for example, is still a hot topic today. Due to its brittleness, cracks and chipping on the edges are generated during cutting (mechanical or laser). Better quality is sought for precision engineering in the microelectronics, display, and watch jewelry industries among others. The mechanical behavior of the material depends on its temperature variation. During heating, thermal expansion occurs and creates tensile stresses in the material. During cooling, the material shrinks and compressive stresses are generated in the material. These two types of thermal stresses coupled with the change of material properties result in residual stresses, which is important to analyze and predict the mechanical strength of the material and "in fine", the product lifetime.

A process combining laser heating and water cooling was even used to cut glass substrate [17]. While scanning, the laser beam is followed by a water jet. A crack is generated at the concordance point between the tensile stresses (laser heating) and compressive stresses (water cooling) acting in opposite directions (**Figure 3**). Laser-material interaction and the crack propagation are analyzed according to thermal stresses generation through numerical simulation in this article [17]. This process results in a clean chipping-free edge.

**Figure 3.** Schematic representation of laser scribe mechanism. (a) Top view: temperature distribution on glass surface. (b) Section view (at A–A): compressive stress induced by heat. (c) Section view (at B–B): tensile stress over compressive stress field. (d) Stress distribution along the central axis at B–B [17].

Although not fully understood, a thermal model was used once again to explain the "stealth dicing" process developed by Hamamatsu [18]. It consists in creating a thin damaged layer inside the material while the surface remains unscathed. It is used for dicing specific micro‐ electronics chips where no debris must be generated and the use of water must be avoided. The authors explain that it is possible to focus the laser beam inside the material when it is partially transparent at the chosen wavelength (1.024 μm wavelength for silicon for instance). At the focal point the intensity is within the range of 108 to 109 W/cm2 . Through heat transfer simulation, it was revealed that the temperature increased and so did the absorption. Me‐ chanical stresses and acoustic waves generated by the temperature raise result in the creation of voids and dislocations in the material structure. The substrate is stuck to a tape. The tape expansion leads to tensile stresses high enough to finally separate the chips.

Other methods are used for more complex processes. In the case of femtosecond lasers, the heating of electrons is decoupled with the lattice. The use of heat transfer equation becomes highly questionable. Different mechanisms can occur depending on the laser intensity. Ionization of the target occurs leaving a mix of charged species repelling each other. Molecular dynamics method is fit to study this kind of phenomenon because it directly studies the movement of atoms and molecules based on Newton equations of motion [19]. However, in this short literature review, we have seen that the thermodynamic and thermo-mechanical phenomena are well suited to describe the dynamics of a laser process. This is also true for ultrashort pulsed laser in some cases. Indeed, depending on the objective, two-temperature model can be implemented for processes where the electron temperature remains low (<104 K). An example of such model is proposed in the following section.

#### **2.3. Application: two-temperature model for picosecond laser formation of LIPPS on copper film**

#### *2.3.1. Context of the study and presentation of the heat transfer model*

The FEM package COMSOL Multiphysics was used to compare the difference between nanosecond and picosecond pulses while modeling laser-copper interaction. The objective here was to qualitatively identify the role of thermal phenomena during the formation of LIPSS (Laser Induced Periodic Surface Structure) on copper by the irradiation of an ultraviolet (266 nm, fourth harmonic) picosecond laser (40 ps). It is an original simulation study that falls within the framework of a former PhD investigation [20]. Indeed, the mechanism of formation of LIPPS is investigated since the 70s and is still under debate. Several mechanisms were proposed and it appeared that the nature of LIPSS and its formation depend on the pulse duration, the wavelength, the material property and state (bulk, film, surface roughness…). While thermal phenomena could be part of the mechanism of LIPSS formation for nanosecond lasers and nonresonant mechanism associated with optical phenomena are responsible for their formation with femtosecond lasers, it is less clear with picosecond lasers [20].

As stated in Section 2.1, the electron heating is decoupled from the lattice heating as soon as the pulse duration is lower than the characteristic time transfer (3 ps for copper). The pulse duration is only one order higher than the electron-lattice time transfer. Therefore, the twotemperature model should be applied in this case (Eqs. (7) and (8)), where electron and lattice are linked together by a coupling parameter *g*(*T*e) (Eq. (10)).

$$C\_\varepsilon(T\_\varepsilon)\frac{\partial T\_\varepsilon}{\partial t} = \nabla(k\_\varepsilon(T\_\varepsilon, T\_\prime)\nabla T\_\varepsilon) - \mathbf{g}(T\_\varepsilon)(T\_\varepsilon - T\_\prime) + \underline{Q} \text{ (heat transfer equation for electrons)}\tag{7}$$

$$\text{R}\_{\text{i}}\text{C}\_{\text{i}}(T\_{\text{i}})\frac{\partial T\_{\text{l}}}{\partial t} = k\_{\text{l}}\nabla^{2}T\_{\text{l}} + \text{g}(T\_{\text{e}})(T\_{\text{e}} - T\_{\text{l}}) \text{ (heat transfer equation for the lattice)}\tag{8}$$

$$k\_{\text{e}} = k\_{\text{e}}^{\cdot} T\_{\text{e}} \; , \; C\_{\text{e}} = C\_{\text{e}}^{\cdot} T\_{\text{e}} \; (\text{for low } T\_{\text{e}}) \tag{9}$$

$$\log(T\_c) = \frac{C\_c}{\tau\_{el}} \text{ (electron-phonon coupling parameter)}\tag{10}$$

$$\mathcal{Q} = \frac{A \times a \times 2 \times F \times \sqrt{\ln(2)}}{\tan \times \sqrt{\pi}} \times \exp\left(-4 \times \ln(2) \times \left(\frac{t - 3tauu}{\tan}\right)^2\right) \ast \exp\left(-a \times x\right) \tag{11}$$

The subscripts "e" and "l" denote, respectively, the electron and the lattice. *T* is the temperature (K), *A* is the optical absorption, α is the linear absorption coefficient (1/m), tau is the pulse duration (s), *F* is the fluence (J/cm2 ), *x* is the depth variable (m), *C* is the thermal heat capacity, *<sup>k</sup>* is the thermal heat conductivity (W/m/K) and *Q* (Eq. (11)) is the heat source (W.m–3). *<sup>C</sup>* ′ *<sup>e</sup>* and *k* ′ *<sup>e</sup>* are, respectively, the derivative of electron heat capacity (mJ/K2 /cm3 ) and the derivative of electron thermal conductivity (W/K2 /m). The material properties used in this investigation are summarized in **Table 1** and were obtained from Hallo et al. [21]:


**Table 1.** Material properties of copper.

At the focal point the intensity is within the range of 108

168 Modeling and Simulation in Engineering Sciences

to 109

simulation, it was revealed that the temperature increased and so did the absorption. Me‐ chanical stresses and acoustic waves generated by the temperature raise result in the creation of voids and dislocations in the material structure. The substrate is stuck to a tape. The tape

Other methods are used for more complex processes. In the case of femtosecond lasers, the heating of electrons is decoupled with the lattice. The use of heat transfer equation becomes highly questionable. Different mechanisms can occur depending on the laser intensity. Ionization of the target occurs leaving a mix of charged species repelling each other. Molecular dynamics method is fit to study this kind of phenomenon because it directly studies the movement of atoms and molecules based on Newton equations of motion [19]. However, in this short literature review, we have seen that the thermodynamic and thermo-mechanical phenomena are well suited to describe the dynamics of a laser process. This is also true for ultrashort pulsed laser in some cases. Indeed, depending on the objective, two-temperature model can be implemented for processes where the electron temperature remains low (<104

**2.3. Application: two-temperature model for picosecond laser formation of LIPPS on copper**

The FEM package COMSOL Multiphysics was used to compare the difference between nanosecond and picosecond pulses while modeling laser-copper interaction. The objective here was to qualitatively identify the role of thermal phenomena during the formation of LIPSS (Laser Induced Periodic Surface Structure) on copper by the irradiation of an ultraviolet (266 nm, fourth harmonic) picosecond laser (40 ps). It is an original simulation study that falls within the framework of a former PhD investigation [20]. Indeed, the mechanism of formation of LIPPS is investigated since the 70s and is still under debate. Several mechanisms were proposed and it appeared that the nature of LIPSS and its formation depend on the pulse duration, the wavelength, the material property and state (bulk, film, surface roughness…). While thermal phenomena could be part of the mechanism of LIPSS formation for nanosecond lasers and nonresonant mechanism associated with optical phenomena are responsible for

their formation with femtosecond lasers, it is less clear with picosecond lasers [20].

As stated in Section 2.1, the electron heating is decoupled from the lattice heating as soon as the pulse duration is lower than the characteristic time transfer (3 ps for copper). The pulse duration is only one order higher than the electron-lattice time transfer. Therefore, the twotemperature model should be applied in this case (Eqs. (7) and (8)), where electron and lattice

ee e e e ( ) ( ( , ) ) ( )( ) (heat transfer equation for electrons) *e e <sup>l</sup> <sup>l</sup>*

¶ (7)

expansion leads to tensile stresses high enough to finally separate the chips.

K). An example of such model is proposed in the following section.

*2.3.1. Context of the study and presentation of the heat transfer model*

are linked together by a coupling parameter *g*(*T*e) (Eq. (10)).

*<sup>T</sup> C T k T T T gT T T Q*

¶ =Ñ Ñ - - +

e

*t*

**film**

W/cm2

. Through heat transfer

Copper films with thicknesses of 0.2, 0.5, and 1 μm were deposited on silicon or glass (not considered in this model) substrate by cathodic magnetron sputtering technique. The laser beam was supposed to have a perfect Gaussian spatial distribution and is represented as heat source *Q*. Thermal losses due to convection or radiation were neglected but can be taken into account as suggested in Ref. [22], especially for longer pulse duration. The optical absorption coefficient and the thermal properties were fixed. Considering the optical penetration depth and the pulse duration, the mesh size of the model and the time step need to be very small. The mesh size was arbitrarily set at 5 nm (optical penetration depth of 12 nm) and the time step was set according to the Neumann criterion (Eq. (12)):

$$
\Delta t < \Delta z^2 \left[ \frac{C(T)}{k(T)} \right]\_{\text{min}} \tag{12}
$$

If the time step is close enough to the minimum requirement, COMSOL will adjust it auto‐ matically according to convergence plots. The mesh size is more important as it is fixed except if adaptive remeshing option is set. It requires a considerable effort for a quad core computer to model this. Fortunately, the spot size is three orders of magnitude bigger than the optical penetration depth. In this case, we can reduce the model to a one-dimensional problem and analyze only the temperature evolution in depth.

#### *2.3.2. Results and discussions*

The single-shot ablation threshold was determined to be around 190 mJ/cm2 by the experiment [20]. The model evaluated the surface peak temperature at 1310 K for 200 nm thick film and 1263 K for 1 μm thick film (**Figure 4**) when a laser fluence of 190 mJ/cm2 is applied. The relaxation is also steeper with a thickness of 1 μm (**Figure 4**) for the pulse duration of 42 ps. The surface peak temperature is higher with a thickness of 0.2 μm and pulse duration of 10 ns. Heat continuity was assumed in this model as there is no information about the interface quality, but it can potentially have an effect on the film temperature especially in the case of nanosecond pulse where the heat penetration depth is around 1 μm. This could cause damage at the copper-substrate interface (delamination). The interface could have a low transmittance and if the substrate is an insulator (glass for instance), the heat could be confined in the layer and the temperature should be higher. This might be the reason why the melting temperature is not reached in the model with a fluence of 190 mJ/cm2 . To be thorough, the interface between the film and the substrate should be modeled as the heat transfer continuity depends on this interface quality. However, with picosecond pulse duration, the heat penetration depth is two orders of magnitude lower than 1 μm but only one order of magnitude lower for 0.2 μm thick film. Thus, the heat transfer continuity might be influent in the case of 0.2 μm thick film but not for 1 μm thick film (**Figure 4**). Also, the material properties were fixed and they are thermodependent in reality, especially when a phase change is involved. Given the values observed for the thermal conductivity, the thermal heat capacity in reference [23], and the fact that the model does not go above the melting point, it was a reasonable assumption, but it might explain the small difference between the model and the experiment. It should be noted that the ablation threshold is evaluated with an error of ±20 mJ/cm2 . The underestimation of the temperature could be due to the absorption as well. A linear absorption coefficient of 0.66 is taken but it could vary because of the ionization of the target [24]. Besides, the properties might be different between the bulk material and thin film. This analysis shows the importance of the initial conditions and simplification hypotheses in a simulation.

and the pulse duration, the mesh size of the model and the time step need to be very small. The mesh size was arbitrarily set at 5 nm (optical penetration depth of 12 nm) and the time

> ( ) <sup>²</sup> ( ) *C T*

*k T* é ù D <D ê ú

If the time step is close enough to the minimum requirement, COMSOL will adjust it auto‐ matically according to convergence plots. The mesh size is more important as it is fixed except if adaptive remeshing option is set. It requires a considerable effort for a quad core computer to model this. Fortunately, the spot size is three orders of magnitude bigger than the optical penetration depth. In this case, we can reduce the model to a one-dimensional problem and

The single-shot ablation threshold was determined to be around 190 mJ/cm2 by the experiment [20]. The model evaluated the surface peak temperature at 1310 K for 200 nm thick film and 1263 K for 1 μm thick film (**Figure 4**) when a laser fluence of 190 mJ/cm2 is applied. The relaxation is also steeper with a thickness of 1 μm (**Figure 4**) for the pulse duration of 42 ps. The surface peak temperature is higher with a thickness of 0.2 μm and pulse duration of 10 ns. Heat continuity was assumed in this model as there is no information about the interface quality, but it can potentially have an effect on the film temperature especially in the case of nanosecond pulse where the heat penetration depth is around 1 μm. This could cause damage at the copper-substrate interface (delamination). The interface could have a low transmittance and if the substrate is an insulator (glass for instance), the heat could be confined in the layer and the temperature should be higher. This might be the reason why the melting temperature

the film and the substrate should be modeled as the heat transfer continuity depends on this interface quality. However, with picosecond pulse duration, the heat penetration depth is two orders of magnitude lower than 1 μm but only one order of magnitude lower for 0.2 μm thick film. Thus, the heat transfer continuity might be influent in the case of 0.2 μm thick film but not for 1 μm thick film (**Figure 4**). Also, the material properties were fixed and they are thermodependent in reality, especially when a phase change is involved. Given the values observed for the thermal conductivity, the thermal heat capacity in reference [23], and the fact that the model does not go above the melting point, it was a reasonable assumption, but it might explain the small difference between the model and the experiment. It should be noted that the ablation

could be due to the absorption as well. A linear absorption coefficient of 0.66 is taken but it could vary because of the ionization of the target [24]. Besides, the properties might be different between the bulk material and thin film. This analysis shows the importance of the initial

*t z*

min

ë û (12)

. To be thorough, the interface between

. The underestimation of the temperature

step was set according to the Neumann criterion (Eq. (12)):

analyze only the temperature evolution in depth.

is not reached in the model with a fluence of 190 mJ/cm2

threshold is evaluated with an error of ±20 mJ/cm2

conditions and simplification hypotheses in a simulation.

*2.3.2. Results and discussions*

170 Modeling and Simulation in Engineering Sciences

**Figure 4.** Lattice surface temperature for different copper thicknesses (0.2 and 1 μm) and different pulse duration (10 ns and 42 ps). The fluence F was set at 190 mJ/cm2 . The time is divided by the respective pulse duration tau in order to have the same scale of observation on the *x*-axis.

The model was tested with pulse duration of 10 ns. It is observed on **Figure 5** that the electron and lattice temperatures are almost identical. It confirms that a one-temperature model is sufficient for nanosecond pulse. For the picosecond model, the electron temperature (3215 K)

**Figure 5.** Electron and lattice surface temperature for pulse duration of 42 ps (top) and 10 ns (bottom) on a 1 μm thick copper target.

is higher than the lattice temperature (1263 K) and there is a delay of the same order as the pulse duration between the two temperatures. It is conformant with the theory which explains that the electrons transfer heat to the lattice during the pulse for the nanosecond regime and it is not the case for the ultrashort regime (especially <500 fs). Here, we have an intermediate regime where the lattice is heated up after the pulse by the electron-lattice heat transfer. The melting point is almost reached. The temperature might be underestimated because of approximations and work hypotheses mentioned earlier. Nevertheless, this model tends to confirm the implication of thermal effects in the formation of LIPPS in the picosecond regime.

This was a single-shot investigation. The damage threshold can be lowered by an accumula‐ tion/incubation effect. It could be interesting to study the thermal implication and the effect of the number of pulses on the delamination of the films. Working under the single-shot ablation threshold with multiple pulses allows a better control of the formation of LIPSS. This is of interest for industrial applications (super hydrophobic surfaces for instance). However, the pulse duration is 9 orders smaller than the pulse repetition rate (40 ps against 10 Hz). The time step could be modulated so that it is short during laser irradiation and longer during relaxation but it would still be difficult to investigate the incubation effect for 1000 pulses at low fluences without the use of a super calculator.

In conclusion, this example shows that COMSOL was able to provide valuable qualitative information about a laser machining mechanism. This model could be improved by lifting the simplifications hypothesis one by one (convection/radiation losses [22], interfaces between solids [25], ionization [26], thermo-dependent properties [27], etc.). Also, with more data (material properties, experimental conditions) it would be possible to have reliable quantita‐ tive information at the expense of more computational power. It should be noted that this method is limited to a nonphase change model and other numerical approaches such as Molecular Dynamic mentioned earlier would be more appropriate otherwise. Overall, Section 2 mentioned how to model situations to better understand the complex physical phenomena involved in a process. The next section presents a more applied approach where direct answers regarding the parameters are desired for process prediction and optimization: the Design Of Experiment.
