Numerical Simulation of Laser Processing Materials: An Engineering Approach

Abstract

The following chapter aims at giving an overview of the use of numerical simulation in the field of laser processing. Indeed, the past two decades saw an increasing demand for lasers in various areas such as healthcare, microelectronics, cartography, optoelectronics, aeronautics, etc. Thus, the comprehension of the laser-material interaction and the removal mechanism became primordial to predict and improve the efficiency of a process. After a nonexhaustive literature review, two simulation approaches (Finite Element and Design Of Experiment, DOE) will be presented to demonstrate the importance of numerical simulation in laser applications.

1. Introduction

As diverse as laser applications are, they have one thing in common: the complexity of the interaction of photons with matter and the multiphysics nature of the phenomena (thermal, fluidic, optical, mechanics, etc.) involved during laser processing whether it is drilling [1], grooving [2], cutting [3], welding [4], estimating ablation threshold limits [5] or simply predicting the thermal effects on the material [6]. Several mechanisms can be involved before, during and after the material ejection and they strongly depend on the laser characteristics (wavelength, pulse duration, beam shape, polarization, etc.) and the process parameters (scanning speed, repetition rate, pulse energy, etc.). The development in computer sciences (calculation capability and software) gave the opportunity to better understand the preponderance or the concomitance of mechanisms depending on the application and the objective of the researcher. Moreover, given the dynamic and short time-scale nature of the laser process, numerical simulation can give an insight of what is happening inside the workpiece, which is delicate by experimental means.

Multiphysics packages such as ANSYS, ABAQUS, ADINA, or COMSOL are effective for the aforementioned issue and for fundamental investigation. In the industry, another semiempirical 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 experimental investigation suffers from constraints such as the availability of the device, the experimenter, the material, the risks, the costs, or the environment.

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.

2. Multiphysics simulation of laser processing

2.1. Introduction

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 of the laser process can be investigated as well.

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 waveguides 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 structuring (LIPSS) on copper film.

2.2. Laser applications and their simulation in the literature

2.2.1. Laser welding

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 absorptivity 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 experimental 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 discontinuities 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 thermo-mechanical 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 conservation is assumed (Eq. (1)).

$\frac{\partial \rho }{\partial t}+\nabla \rho {\displaystyle \overrightarrow{u}}=0\text{ (mass conservation equation)}$E1

$\rho \frac{\partial {\displaystyle \overrightarrow{u}}}{\partial t}+\rho {\displaystyle \overrightarrow{u}}\nabla {\displaystyle \overrightarrow{u}}=-\rho \nabla \text{P}+\eta \Delta {\displaystyle \overrightarrow{\text{u}}}\text{ (Navier-Stokes)}$E2

$P_0+P_{\text{v}}+P_{\text{L}}=2{\rm K}\sigma$E3

where ρ is the material density, t is the time variable, u is the velocity field, P₀ 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:

$\frac{\partial \rho H}{\partial t}-\nabla k\nabla T+\nabla \rho {\displaystyle \overrightarrow{u}}H=Q$E4

$H=L+C_{\text{p}}\times T$E5

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 modification 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 10⁸ W/cm².

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 CO₂ 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 modification over the phase change temperature range (Eq. (6)).

$h(T)=\text{flc2hs}(T-T1,T2-T1)$E6

where T1 is the temperature where the phase change begins, T2 is the temperature where it ends, and flc2hs is a smoothed Heaviside function used in COMSOL. T1, T2, 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

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 microelectronics 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 10⁸ to 10⁹ W/cm². Through heat transfer simulation, it was revealed that the temperature increased and so did the absorption. Mechanical 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 (<10⁴ 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 two-temperature 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_e(T_e)\frac{\partial T_{\text{e}}}{\partial t}=\nabla (k_{\text{e}}(T_{\text{e}},T_l)\nabla T_{\text{e}})-g(T_{\text{e}})(T_{\text{e}}-T_l)+Q\text{ (heat transfer equation for electrons)}$E7

$C_{\text{l}}(T_{\text{l}})\frac{\partial T_{\text{l}}}{\partial t}=k_{\text{l}}{\nabla }^2{T}_{\text{l}}+g(T_{\text{e}})(T_{\text{e}}-T_{\text{l}})\text{ (heat transfer equation for the lattice)}$E8

$k_{\text{e}}=k_{{}_{\text{e}}}^{\text{'}}{T}_{\text{e}},C_{\text{e}}=C_{\text{e}}^{\text{'}}{T}_{\text{e}}\text{ (for low }{T}_{\text{e}}\text{)}$E9

$g\text{(}{T}_{\text{e}}\text{)}=\frac{C_{\text{e}}}{{\tau }_{\text{el}}}\text{ (electron-phonon coupling parameter)}$E10

$Q=\frac{A\times \alpha \times 2\times F\times \sqrt{\ln (2)}}{\text{tau}\times \sqrt{\pi }}\times \exp \left(-4\times \ln (2)\times {\left(\frac{t-3\text{tau}}{\text{tau}}\right)}^2\right)*\exp \left(-\alpha \times x\right)$E11

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/cm²), x is the depth variable (m), C is the thermal heat capacity, k is the thermal heat conductivity (W/m/K) and Q (Eq. (11)) is the heat source (W.m^–3). ^{${C^{\prime }}_e$} and ${k^{\prime }}_e$ are, respectively, the derivative of electron heat capacity (mJ/K²/cm³) and the derivative of electron thermal conductivity (W/K²/m). The material properties used in this investigation are summarized in Table 1 and were obtained from Hallo et al. [21]:

Density (kg.m^–3)	8960
Number of electrons per atom Z	1
Optical penetration depth ls (m)	12.5 × 10^–9
lattice heat capacity Cl (J/K/m^3)	395
Derivative of electron heat capacity Ce’ (J/K^2/m^3)	70
Lattice thermal heat conductivity kl (W/m/K)	390
Derivative of electron heat conductivity ke’ (W/K^2/m)	1.34
Electron-lattice exchange rate g (W/K/m^3)	34 × 10^15
Melting temperature (K)	1360

Table 1.

Material properties of copper.

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²{\left[\frac{C(T)}{k(T)}\right]}_{\min }$E12

If the time step is close enough to the minimum requirement, COMSOL will adjust it automatically 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/cm² 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/cm² 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/cm². 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 thermo-dependent 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/cm². 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.

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/cm². 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) 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.

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.

This was a single-shot investigation. The damage threshold can be lowered by an accumulation/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 quantitative 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.

3. The Design Of Experiment methodology

3.1. Introduction

The Design Of Experiment methodology is based on statistical analysis and provides a semiempirical model in a similar way than limited development. Indeed, a response (laser hole drilling depth, for instance) can be described by a polynomial function of several parameters in a restricted domain through multilinear regression. Opposed to the Change One Separate factor at a Time (COST) method, DOE saves time and resources (human, machine, costs, sample, etc.) and is more relevant when parameters have a coupled effect on the response. It is a powerful tool to rule out noninfluent parameters (screening), describe, test the robustness or even optimize a process with a minimum of experimental trials. It is widely used in the industry to develop new products/process or improve existing ones, improve the quality, reduce the production costs, or the impact on the environment and even evaluate the influence of perturbations on a process (temperature, humidity, etc.) [28, 29].

The following sections propose to explore the use of this method in laser-based applications and to discuss the study of dicing silicon carbide substrates in the microelectronics industry.

3.2. The use of Design Of Experiment in the literature

The following paragraph aims at explaining the key steps and objectives of the Design Of Experiment methodology and cannot be taken as a full course on the subject. A literature review reveals that many authors have used Taguchi tables or the more conventional Response Surface Methodology (RSM). It is important to note that several strategies exist and they must be used depending on the objective of the researcher and the means/costs involved. The choice of strategy is strongly dependent on the researcher’s experience and knowledge, especially when it comes to defining the control parameters, identifying the noise parameters (the one we are submitted to), the response he wants to characterize/measure, and the range of variation of the parameters (domain of study). Once these steps have been completed, whatever the employed strategy, it needs to be addressed with great care and in a sequential manner.

3.2.1. The screening

The parameters are coded between –1 (minimum value) and +1 (maximum value) so that the scale of variation becomes the same between the parameters and it is easier to make a hierarchy of their degree of influence. In RSM, the first step is the analysis of the center of the domain. All parameters are set at 0 (average value) in their coded form. The experiment is replicated three times. If the value is two orders of magnitude larger than the standard error, it is possible to carry on the study of the response. It means that the response obtained is not due to perturbation/noise. Otherwise, there might be a problem with the measurement tools and most likely that the domain of study is not large enough. It must be noted that the number of replication of the center point depends on the number of parameters and levels and needs to be increased when a quadratic model is envisaged.

The following step, the screening, is crucial and can already provide significant information such as: determining which parameters are influent and create a hierarchy, postulating a linear model without interaction, providing a predictive model (qualitatively), and determining if the domain was established properly. In some cases, it is already noticeable that interactions between parameters are influent and that nonlinearity is expected. Remember, if the objective was to identify among a lot of parameters which one were the most significant, then there is no need to carry on. A screening postulates a linear model without interaction.

A proper way to conduct the screening is to use the Hadamard matrices. These are square matrices (H2, H4, H8, H16, etc.) and correspond to a system of equations. With three parameters for instance, H4 is recommended as three equations are required to solve the three unknown parameters and one degree of freedom is available (for the average value). With four parameters, H8 is recommended because H5 does not exist. The four remaining trials are not lost since Hadamard matrices are contained into full factorial or composite designs which are used to obtain, respectively, a linear model with interactions and a quadratic model.

3.2.2. Prediction and optimization

For prediction or optimization objectives, a full factorial design (linear model + interaction) or a composite design (quadratic or cubic model) is carried on to have a relevant predictive model that will allow setting the parameters values according to the desired response. When a nonlinear model is involved, an optimization step is possible. It is even more useful when several responses are studied and that a compromise must be found. A function called “desirability function” is used especially to optimize multiple responses simultaneously according to the target of the scientist. Each response Y_i is transformed into a dimensionless factor bounded by 0>Di<1 where a higher “Di” indicates a more desirable response value. Three cases arise. “Nominal-the best” case is when the scientist desires to reach a certain response value. “Larger-the best” is when the scientist wants to maximize a response value (laser welding depth of penetration for instance at the highest possible speed) and the “smaller-the better” case is when he wants to minimize a response (HAZ for instance). That desirability function is usually treated with statistical software according to the scientist’s objective (minimize, maximize, targeted value). Ref. [30] describes a good example of prediction and optimization of laser butt welding. The bead width, depth of penetration, and the tensile strength of the joint are analyzed. The desirability function is combined with the use of Taguchi tables in order to optimize the parameters (beam power, travel speed, and focal position) and improve the quality of the weld. In this chapter, the quality of the weld is optimal when the bead width is minimized and the tensile strength and depth of penetration are maximized. Then, the desirability function is set accordingly. For this configuration, the laser beam power was found to be the most significant parameter. The interaction effects between the beam power and travel speed, the travel speed, and the focal position were found to be significant as well when using helium shielding gas (Figure 6). The interaction effect between beam power and travel speed was found to be most influent followed by the beam power, the travel speed and the focal position when other gases were used.

Figure 6.

Factor effects on composite desirability values for helium gas [30].

Figure 7.

Desirability function analysis and results obtained with optimal parameters. (a) Surface view of laser scribed sapphire, (b) wall section view [31].

In another example, the authors in Ref. [31] wanted to dice sapphire substrates used in the fabrication of microelectronics chips with a Nd:YAG nanosecond laser (second harmonic 532 nm). They aimed at a depth equivalent to 1/3–1/4 of the wafer thickness (108–148 μm in this case) to be able to subsequently break it and separate the dies. In the meantime, they must choose carefully the laser parameters in order to minimize the groove width because of economical and throughput reason. The studied parameters were the pulse energy (A), the scanning velocity (B), and the number of passes (C). It was found that A, B, and C were influent on the depth as well as the interaction terms AB, AC, BC, and the quadratic terms A² and B². The effects of pulse energy and scanning times have slight and nearly equal effects while speed was the dominant parameter especially when the laser is scanned several times with high pulse energy. For the width, the dominant parameter was the pulse energy. A combination of low scanning velocity (0.5–0.59 mm/s), low pulse energy (150 μJ/pulse), and multiple passes (3) was found based on desirability analysis (Figure 7) to obtain deep (148 μm) and narrow grooves (19 μm) (Figure 7a and b).

In this example [31] and others [12, 32], Box Benhken Designs (BBD) were used. They are an alternative to composite designs because they contain fewer experiments than the latter. However, there might be a loss of accuracy in the final model; they do not contain the screening and a nonlinear model is postulated from the start. In the chapters mentioned above, it is suggested that parametric studies were performed before doing the BBD. Then the authors knew the model would be nonlinear. Overall, the risk of using BBD must be assessed carefully.

A problem of laser welding of AISI304 stainless steel is presented in [12]. The depth of penetration of the weld and the bead width characterize the geometry and are related to the mechanical quality of the weld. These two responses are studied according to the beam power, the welding speed, and the beam angle using a BBD. It was found that the beam power was the most significant parameter followed by the welding speed. This is conformant with Ref. [30] mentioned earlier. With an optimal combination of parameters 1250 W, 750 mm/min, and a beam angle of 90°, depth of penetration of 1.5 mm and beam width of 1.3 mm were achieved. The results were compared to finite element simulation and were discussed in Section 2.2.

In Ref. [32], BBD was used to predict the geometry (depth and width) of CO₂ laser-machined micro-channel in glass. Three different models using Artificial Neural Network (ANN) were implemented through Labview^®. It was found that one of them gave a better prediction precision than DOE and the two others had a greater average percentage error. Overall, the ANN strategy was deemed useful for prediction and laser parameters optimization. It was found that the laser beam power had a positive effect on the channel width and depth increase. The pulse repetition rate had a negative effect on both responses. The traverse speed had a negligible effect on the channel width and tended to decrease the channel depth when increasing.

3.2.3. Prediction and robustness

The study of the robustness is another objective and is dear to the industrials for better process control and repeatability. It can be performed through the analysis of the error (for each experiment) as a response while carrying on RSM. Taguchi method is a particular case of the Design Of Experiment and is relevant for robustness study. The tables of Taguchi are mostly fractional designs. In Taguchi methodology, the signal to noise (S/N) ratio is analyzed according to the following equation:

$\text{S/N}i=-10\log \left[\frac{1}{n}{\displaystyle {\displaystyle \sum _{i=1}^n}{y}_i^2}\right]$E13

$\text{S/N}i=-10\log \left[\frac{1}{n}{\displaystyle {\displaystyle \sum _{i=1}^n}\frac{1}{y_i^2}}\right]$E14

$\text{S/N}=-10\log \left[\frac{1}{n}{{\displaystyle {\displaystyle \sum _{i=1}^n}(y_i-m)}}^2\right]$E15

where n is the number of trials for each experiment, y_i is the response of the ith experiment, and m is the target the scientist wants to reach. It allows control and the reduction of variability of the response. The “smaller-the better-characteristic” (Eq. (13)) is used when the scientist wants to minimize a response as opposed to the “larger-the better characteristic” (Eq. (14)). The “nominal-the-best” characteristic (Eq. (15)) is about reaching a specified target. It is not the same as the desirability function (Section 3.2.2) since this is an analysis of the signal-to-noise ratio and this is what makes the Taguchi method more relevant for robustness study.

The methodology is a bit different than RSM [30]: (1) identify the control factors and the interactions, (2) identify the levels of each factors, (3) select the appropriate orthogonal array (OA), (4) assign the factors/interactions to the columns of the OA, (5) conduct the experiments, (6) analyze the data and determine the optimal levels, (7) perform the confirmation experiments. The two first steps are already tricky and imply that the researcher must have a prerequisite knowledge or have performed earlier experiments to guide him. For steps 3 and 4, the researcher has to choose the orthogonal array and how to fill the columns, thanks to steps 1 and 2 and thanks to the “interaction graphs.” The Taguchi tables are indeed limited to a certain type of model with a restricted number of parameters and levels. They have to be filled according to the “interaction graphs.” They show the main factors and the interaction that are possible to study for each Taguchi table. These graphs contain other criteria as well: they show how difficult it is to change the level of parameters [33]. For processes where the temperature is a key factor, it takes more time to cool down an oven than to heat it for instance. An experiment in this case can have an effect on the following one. Therefore, the sequence of experiment must be chosen carefully. Taguchi method has been successfully used in numerous studies [34, 35].

After each batch of experiments, the Analysis Of Variance (ANOVA), the analysis of multilinear regression coefficient and the study of residuals must be carried on. These are statistical means to validate the models. Once the final model is found to be statistically adequate, new experiments (not included in the original design) must be executed and compared to the model for its final validation. Another objective of the work cited as Ref. [32] was to compare DOE and ANN methods and show the interest in developing new modeling techniques with the same robustness than DOE while being more easily or successfully applied. This should be kept in mind as the way to go for future simulation method development. The overall Design Of Experiment methodology is summarized by Figure 8.

Figure 8.

Response surface methodology diagram. Yellow designates the analysis and decisions depending on the researcher’s experience, squares designate the experimental work, and the red one designates a turning point where lot of information is gathered and will decide on the next steps.

3.3. Application: optimization of laser parameters for dicing silicon carbide substrate

3.3.1. RSM study

The following study is extracted from Ref. [36]. In the microelectronics industry, laser is used to replace the conventional method of blade dicing for hard and chemically inert substrates. Laser is used to scribe and then break the substrate in order to separate electronic chips. A Diode Pumped Solid State (DPSS) Nd:YAG laser operating at a wavelength of 355 nm (3rd harmonic), a repetition rate of 40 kHz with a pulse duration around 90 ns, and a spot size around 5 μm is used. The laser source is embedded in an enclosed station as it is destined for production and is operating in white room environment. The substrate of interest is a three inches diameter silicon carbide wafer and is 360 μm thick. The scribe and break method requires a scribing depth of one forth to one third of the wafer thickness. However, the breaking step becomes difficult for hard substrates such as SiC and because the breaking efficiency depends on the leverage (dies size). The groove needs to go deeper and the objective is to optimize this depth according to the laser parameters: pulse energy (A), number of passes (B), defocus (C), and scanning speed (D). As seen in Section 3.2, other parameters are influent on a laser process but these are the main parameters that an operator can customize on this station. The parameters and their values are reported in Table 2.

Parameter (Pi)	Identifier	Level –1 (low)	Level 1 (high)	Level 0 (average)	∆Pi
Energy (μJ/pulse)	A	30	140	85	110
Passes	B	1	5	3	4
Defocus (μm)	C	–15	–65	–40	–50
Speed (mm/s)	D	10	48	29	38

Table 2.

Parameter values, identifiers, and study domain.

Following the procedure illustrated by Figure 8, the experiment performed at the center of the domain is replicated seven times (in italic in Table 3) and reveals that the response can be studied because the mean response value (157.9 μm) is two orders greater than the standard deviation (1.97). In other words, the response obtained is due to the parameters and not due to noise factors.

Sequence	Energy (A)	Passes (B)	Defocus (C)	Speed (D)	Depth (Y)
1	–1	–1	–1	–1	73.8
2	1	–1	–1	–1	210.8
3	–1	1	–1	–1	101
4	1	1	–1	–1	286.2
5	–1	–1	1	–1	72
6	1	–1	1	–1	181.2
7	–1	1	1	–1	100
8	1	1	1	–1	261
9	–1	–1	–1	1	30
10	1	–1	–1	1	127.2
11	–1	1	–1	1	65
12	1	1	–1	1	180
13	–1	–1	1	1	47.2
14	1	–1	1	1	122
15	–1	1	1	1	68.4
16	1	1	1	1	165.8
17	–α	0	0	0	85.6
18	A	0	0	0	189.6
19	0	–α	0	0	124.4
20	0	α	0	0	168.6
21	0	0	–α	0	156.4
22	0	0	α	0	143.2
23	0	0	0	–α	201.6
24	0	0	0	α	123.2
25	0	0	0	0	158.8
26	0	0	0	0	158.6
27	0	0	0	0	158.4
28	0	0	0	0	159
29	0	0	0	0	156.2
30	0	0	0	0	154.4
31	0	0	0	0	160.2

Table 3.

Central composite face-centered in a cube design (α=1) and measured response values.

Source	Sum of squares	Degree of freedom	Mean square	F value	p-Value
Model	40252.44	6	6708.74	2024.44	<0.0001
A	28752.02	1	28752.02	8676.25	<0.0001
B	2298.42	1	2298.42	693.57	<0.0001
C	220.5	1	220.05	66.54	<0.0001
D	7176.02	1	7176.02	2165.45	<0.0001
AD	1740.5	1	1740.5	525.22	<0.0001
BD	64.98	1	64.98	19.61	<0.0031
Curvature	2660.03	1	2660.03	802.69	<0.0001
Residual	23.20	7	3.31
Lack of fit	0.02	1	0.02	5.178E-3	0.9450
Pure error	23.18	6	3.86
Cor total	42935.67	14

Table 4.

ANOVA table for the screening (Hadamard design).

Table 3 shows the full design. The screening is set as a Hadamard H8 (because four factors are involved and H5 does not exist) and is shown in normal font in Table 3. The ANOVA results (Table 4) are used to determine the relevancy of a linear model. We focus on the regression coefficients and the p-value associated to the lack-of-fit and each parameter. The p-value (Fisher test) must be lower than 0.05 for the item to be declared significant. Thus, in this case, the lack of fit is not significant. The multilinear regression coefficients are all low and the difference of mean values at the central point for the model and the experiment is greater than the standard deviation. A linear model is obviously not adequate. A linear model taking into account the interactions between parameters will not be adequate either. Indeed, the curvature is deemed significant by the ANOVA. It means that a nonlinear model will be adequate.

Source	Sum of squares	Degree of freedom	Mean square	F value	p-Value
Model	1.02E+5	16	6372.54	624.62	<0.0001
A	5408.0	1	5408.0	530.07	<0.0001
B	9220.82	1	9220.82	903.79	<0.0001
C	269.12	1	269.12	26.38	0.0002
D	3073.28	1	3073.28	301.23	<0.0001
AB	1232.01	1	1232.01	120.76	<0.0001
AC	529.0	1	529.0	51.85	<0.0001
AD	2704.0	1	2704.0	265.04	<0.0001
BD	207.36	1	207.36	20.32	0.0005
CD	216.09	1	216.09	21.18	0.0004
A^2	845.87	1	845.87	82.91	<0.0001
B^2	217.46	1	217.46	21.31	0.0004
C^2	88.93	1	88.93	8.72	0.0105
D^2	118.11	1	118.11	11.58	0.0043
ABD	222.01	1	222.01	21.76	0.0004
A^2D	149.65	1	149.65	14.67	0.0018
AB^2	145.6	1	145.6	14.27	0.002
Residual	142.83	14	10.2
Lack of fit	119.66	8	14.96	3.87	0.0581
Pure error	23.18	6	3.86
Cor total		30

Table 5.

ANOVA table for the cubic model.

A central composite face-centered in a cube design was chosen (Table 3) to model a quadratic, possibly a reduced cubic model. A third parameter α is introduced in a central composite design. Its value depends on the number of parameters and levels and is chosen to cope with the statistical degradation of the accuracy of the model. Indeed, nonlinear models require the use of three levels which make 81 experiments to execute for four parameters (3⁴). The use of this third value α allows decreasing the number of experiments compared to a full design while keeping a certain statistical accuracy (although not as good as a full design).

Finally, the ANOVA results (Table 5) reveal that the model is significant and the lack-of-fit is not (p-value close to 0.05). All coefficients shown in the ANOVA table are significant; the others are not, hence the term of “reduced” cubic model. The multilinear regression coefficients are close to 1 (R² = 0.9986, R² adjusted = 0.997, R² predicted = 0.9805). The reduced cubic model is thus significant and the final mathematical equation (Eq. (16)) in terms of coded factors is given using the least-squares method.

Y = 156.89 + 52A + 22.63B – 3.87C -– 39.2D + 8.78AB – 5.75AC – 13AD – 3.6BD + 3.68CD – 18.05A² – 9.15B² – 5.85C² + 6.75D² – 3.73ABD + 9.17A²D + 9.05AB²

E16

The model is determined to be statistically adequate. The study of the residuals did not reveal major inconstancies. It is necessary to check if the model is scientifically adequate and if it can predict response values with additional experiments.

3.3.2. Experimental validation

From the response equation (Eq. (16)), it is clear that the pulse energy (A) and the scanning speed (D) are the most influent parameters on the scribing depth. Heating a target via laser irradiation depends on the amount of energy brought to the target and the interaction duration. This is why they are the most dominant factors and why the interaction term AD appears in the equation. Square terms coefficient A² and D² are conformant with theory as well. Indeed, the depth increases with increasing energy input until saturation. It is the classical logarithm trend observed in the literature [37]. The scribing depth also decreases with increasing speed until saturation. Interaction terms such as AC (evaluation of the fluence, Eqs. (18) and (19)), AD (effect of fluence over time), and CD (laser/material interaction over time) appear naturally in the response equation because they are linked by Eqs. (17)–(19):

$F=\frac{E}{\pi \times \omega (z)²}=\frac{A}{\pi \times \omega (C)²}$E17

$O=1-\frac{v}{\omega (z)\times f}=1-\frac{D}{\omega (C)\times f}$E18

$\omega (z)={\omega }_0\sqrt{\left(1+{\left(\frac{M²\lambda z}{\pi {\omega }_0^2}\right)}^2\right)}={\omega }_0\sqrt{\left(1+{\left(\frac{M²\lambda C}{\pi {\omega }_0^2}\right)}^2\right)}$E19

where z is the defocus (m), O is the overlap ratio, f is the repetition rate (Hz), v is the scanning speed (mm/s), λ is the wavelength (nm) and ω₀ (cm) is the spot radius at the focal point. The effect of the number of passes B varied depending on the delay between two passes and thus on the scanning speed. The effect saturated around three passes for pulse energies below 85 μJ and tended to saturate around five passes above this energy threshold. The number of scans effect saturates as the aspect ratio (depth/width) increases because there is less margin for debris ejection. Debris disturbs the beam propagation at the bottom of the trench and successive passes become less effective. This inconvenience can be balanced depending on the pulse energy and the scanning speed because they rule the laser-material interaction efficiency. This is why interaction terms between the number of passes (B), the pulse energy (A), and the scanning speed (D) appear in the response equation.

Additional experiments were performed and compared to the predicted values given by the model. A reasonable correspondence was found between predicted and experimental values as observed in Figure 9.

Figure 9.

Graph showing the variation of depth versus the energy at different scanning speeds. The frequency was fixed at 40 kHz.

This model was successfully used to determine optimal parameters for an objective of 180 μm deep scribe in silicon carbide substrate.

4. Conclusion

Whether it is Finite Element Methods, Monte Carlo or Molecular Dynamics, numerical simulation gives the opportunity to better comprehend the physics of laser-material interaction and in many cases to predict results and optimize the parameters for laser processing. In this chapter, numerical simulation was proven to be valuable to both fundamental science and laser applications in the industry. Computer calculation improvement makes it possible to model complex problems when multiphysical phenomena are coupled and allowed important achievements in the field of laser processing. Simplification assumptions are still required but with time and research efforts, models will become more and more realistic giving the possibility of getting a snapshot inside complex physical mechanisms that cannot be observed through experimental means.