Nanosecond Laser Ablation: Mathematical Models, Computational Algorithms, Modeling

Abstract

The basic mathematical models, computational algorithms, and results of mathematical modeling of various modes of laser action on metals are considered. It is shown that for mathematical description and analysis of the processes of laser heating, melting, and evaporation of condensed media, various theoretical approaches are used: continuum, kinetic, atomistic, etc. Each of them has its own field of applicability, its advantages, and disadvantages. Mathematical description of ns-laser ablation is usually carried out within the framework of continuum approach in the form of hydrodynamic models that take into account reaction of irradiated material to varying density, pressure, and energy both in the target and in the vapor-gas medium. Within the framework of continuum approach, a multiphase, multifront hydrodynamic model and computational algorithm were constructed that were designed for modeling ns-PLA of metal targets embedded in gaseous media. It is shown that proposed model and computational algorithm allow to carry out the simulation of interrelated mechanisms of heterogeneous and homogeneous evaporation of metals manifested as a series of explosive boiling. Modeling has shown that explosive boiling in metals occurs due to the presence of a near-surface temperature maximum. It has been established that in ns-PLA, exposure regimes can be realized in which a phase explosion is the main mechanism of material removal. The verification of reliability of obtained results was carried out by comparing experimental data and calculations with atomistic models.

1. Introduction

The pulsed laser ablation (PLA) of condensed media has been intensively studied over the past few decades [1, 2]. The increased interest in PLA is determined by the increasing possibilities of its use in a variety of applications, beginning with the already traditional ones: microprocessing [3, 4], pulsed laser deposition (PLD) [5, 6], laser-induced breakdown spectroscopy (LIBS) [7, 8], and new rapidly developing areas: production of nanomaterials [9, 10], surface nanostructuring [11], chemical, and physical synthesis [12]. The total effect of nanotechnology and photonics led to the emergence of a new direction—laser synthesis of colloids [13], which draws general attention to its extensive applicability, primarily in biomedicine [14, 15].

Numerous applications make PLA an attractive direction for fundamental investigations. Despite extensive studies of fundamental properties of laser ablation performed earlier, a number of important physical phenomena still remain insufficiently well studied and understood. Previous studies have established that the nature of interaction of laser radiation depends both on the modes of action: wavelength [16, 17], duration [18, 19], and laser pulse intensity [20, 21] and on thermophysical and optical properties of the target [22, 23], presence of the surrounding gas [24] and its pressure [25, 26]. The greatest differences in the physical mechanisms of laser ablation of metals are observed between short (ns) and ultrashort (ps, fs) pulsed modes [18, 19].

In ultrashort range (fs, pc) of influence, laser radiation freely reaches surface of the target. Absorption of laser radiation by a degenerate electron gas followed by a slowed-down energy exchange between electron and phonon components leads to strong deviation from locally thermodynamic equilibrium of the system as a whole. As a result, laser ablation and its accompanying processes develop after the end of the pulse.

Laser ablation in the nanosecond range is a more complex phenomenon involving many interrelated processes both during action and after the end of laser pulse. Such processes include laser target heating, heterogeneous, and homogeneous phase transformations, taking place in the evaporated matter, formation and expansion of plasma plume, heat transfer, laser and intrinsic plasma radiation transfer, generation and propagation of shock waves, and contact boundaries in the vaporized matter and the surrounding gas environment. In contrast to ultrashort regime, in nanosecond range, with a certain choice of parameters of the action, two experimentally observed and explored phenomena arise—volumetric boiling (phase explosion) of liquid phase of the target and formation of laser plasma in vaporized matter and surrounding gas.

Putting an irradiated target into an external gas environment, which is typical for most PLA applications, significantly complicates the situation. In these cases, the long-lived processes of laser-plasma plume and associated generation and propagation of interacting fronts of shock waves and contact boundaries in the vapor-gas medium play an important role in overall picture of laser ablation [27]. The presence of such a large number of interconnected physical processes creates additional difficulties in determination and investigation of the basic mechanisms of ablation. The information obtained by the methods of instrumental diagnostics is not sufficient due to the lack of data on fundamental phenomena associated with rapidly changing energy (thermal and laser radiation), hydrodynamic fields, and the kinetics of heterogeneous and homogeneous phase transformations in the solid and liquid phases, ionization of vaporized matter, and gas environment. At the same time, the understanding of fundamental physics of internal structure of plasma plumes and their spatiotemporal evolution in the process of plasma expansion in the background gas at atmospheric pressure is of decisive importance for many engineering applications. For this reason, nanosecond laser ablation continues to be an area of active research in which mathematical modeling plays an increasing role [28–31].

Any simulation begins with a choice of mathematical model, construction and development of which is given a paramount importance in the computing experiments. For theoretical description and analysis of PLA process of condensed media, various theoretical approaches are used: continuum, kinetic, atomistic (molecular dynamics, etc.). Each of them has its own field of applicability, its advantages, and disadvantages.

Atomistic models allow us to conduct research at the atomic level and obtain fundamental knowledge about structure, thermodynamic, and mechanical properties of crystalline materials [32, 33], about physical mechanisms of various processes [34, 35], including the kinetics of heterogeneous and homogeneous phase transitions [36, 37]. The basic methods of atomistic modeling—molecular dynamics (MD) and Monte Carlo (MC), which use as a rule, semiempirical interaction potentials, operate with tens and hundreds of millions of atoms, and allow calculations in the time range of nanosecond duration.

However, even with the use of parallel computer platforms, computational costs are enormous, and the space-time scales inherent for PLA processes are beyond the limits of accessibility for atomistic modeling methods. Therefore, in spite of constant progress in the field of designing interatomic potentials and increasing the power of computing systems, the final overcoming of computational constraints is hardly achievable, and continuum models will always remain relevant.

Continuous models based on the equations of continuous medium mechanics are realized, as a rule, in the form of hydrodynamic models [18, 24, 28, 30, 31, 38–42] and use the minimum of information and operate with average values of physical characteristics calculated on infinitesimal volume. The methods for solving them are more compact, they have higher accuracy and a relatively small amount of computation. The main shortcomings of continuum approach are manifested in the absence of the possibility of direct investigation of elementary processes in materials and limited possibilities of mathematical description of homogeneous mechanisms of phase transitions of the first kind and calculation of thermophysical, thermodynamic, optical, and other characteristics of matter in a wide range of parameters. These problems are much easier and more fully solved within the framework of atomistic modeling, the results of which can be used as input parameters in meso- and macrolevel models.

In this paper, the application of continuum approach to modeling in preplasma regime of processes dynamics and of main mechanisms of ns-PLA of metal target (Al) in air is considered. The mechanisms of heterogeneous and homogeneous phase transitions interacting with each other are analyzed in detail. Primary attention is paid to model the dynamics of phase explosion of liquid phase of aluminum and the expansion of its fragments in the air, since explosive boiling is considered to be one of the most effective thermal mechanisms of ns laser ablation of materials. Various aspects of this problem have been studied in a number of theoretical and experimental studies [43–56], but there is still no consensus on the mechanism of the phase explosion in metals. In order to obtain detailed information on interaction of heterogeneous and homogeneous mechanisms and data on laser plume morphology, simulation of laser heating, melting, surface evaporation, and evolution of plume in the vapor-gas medium is performed within the framework of new hydrodynamic model with temperature dependences of material properties of the target and explicit tracking of interphase boundary fronts, contact boundary, and shock wave. The release of liquid phase fragments into the atmosphere as a result of phase explosion is modeled by the procedure for introducing a quasinucleus of a new phase (vapor) in the region of near-surface maximum of temperature reaching the value T_max  ≈ 0.9T_cr at time of reaching the maximum permissible overheating of metastable liquid phase.

2. Hydrodynamic model

Laser ablation of Al target placed in the air atmosphere is considered. In preplasma regimes of the action in ablation, the processes in three phase states (solid, liquid, and vapor) and two states (perturbed, unperturbed) of external gaseous medium are taken into account. Laser radiation propagates from right to left. The air for the selected action mode is completely transparent to the laser flow, which is partially reflected from the metal surface and partially absorbed in the near-surface layer of the target.

The mathematical description of ns laser ablation within the continuum approach is realized in the form of hydrodynamic models that allow one to take into account the reaction of condensed (target) and vapor-gas (vaporized matter, gas) media to varying density, pressure, and energy. The processes in each medium are described by a system of nonstationary equations of gas-hydrodynamics supplemented by equations of energy with thermal conductivity, equation of laser radiation transfer, and corresponding equations of state. The 1-D approximation is used for spatial variables. For laser action regimes under consideration, conditions of locally thermodynamic equilibrium (LTE) are assumed for all processes both in the irradiated target and in the vapor-gas medium. Accordingly, hydrodynamic model is formulated in one-temperature approximation.

Phase states and vapor-gas medium are separated among themselves by moving interphase boundaries solid-liquid Г_sℓ (t), liquid-vapor Г_ℓv (t), and by the fronts of contact boundary Г_vg (t) and shock-wave Г_sh,g (t). The right-hand boundary Г_g (t), running along unperturbed gas, is declared moving in order to improve economic efficiency of the computational algorithm. Schematically, their position and direction of motion are shown in Figure 1 .

Figure 1.

Scheme of spatial position of the phases and the direction of motion of interphase, contact boundaries, and the front of shock wave.

In order to obtain complete information on the kinetics of phase transformations, the morphology and dynamics of laser plume, all moving fronts and boundaries in the course of the solution are subjected to explicit tracking using appropriate relationships, which are simultaneously the boundary conditions for hydrodynamic and energy equations.

Complete system of equations is represented as:

${\left[\begin{array}{l}\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u\right)}{\partial x}=0,\\ \frac{\partial \left(\rho u\right)}{\partial t}+\frac{\partial \left(\rho u^2\right)}{\partial x}+\frac{\partial P}{\partial x}=0,\\ \frac{\partial \left(\mathit{\rho \epsilon }\right)}{\partial t}+\frac{\partial \left(\rho \mathit{u\epsilon }\right)}{\partial x}=-\left(P\frac{\partial u}{\partial x}+\frac{\partial W_T}{\partial x}+\frac{\partial G}{\partial x}\right),\\ \frac{\partial G}{\partial x}+{\kappa }_L\left(\rho ,T\right)G=0,\\ W_T=-\lambda \left(T\right)\frac{\partial T}{\partial x}\\ P=P\left(\rho ,T\right),\epsilon =\epsilon \left(\rho ,T\right)\end{array}\right]}_{k}k=s,\ell ,\upsilon ,g$ E1

$t>0,{\Gamma }_s<x<{\Gamma }_{\mathit{s\ell }}\left(t\right)<x<{\Gamma }_{\mathit{\ell v}}\left(t\right)<{\Gamma }_{\mathit{\upsilon g}}\left(t\right)<{\Gamma }_g\left(t\right)$

Here, ρ, u, ε, T, and P are the density, gas dynamic velocity, internal energy, temperature, and pressure of the substance, respectively, κ_L and G are the absorption coefficient and laser radiation fluence, W_T is heat flux density, and λ is coefficient of thermal conductivity. Indices s, ℓ, v, g denote the belonging of the quantities to solid, liquid, vapor, and air, respectively. In the condensed phase, the value ε_k has the meaning of enthalpy of liquid and solid phases H_k .

2.1. Boundary conditions

In the hydrodynamic models that easily combine with the kinetic ones, the heterogeneous mechanisms of phase transitions: melting/crystallization and evaporation/condensation are described naturally. A distinctive feature of heterogeneous phase transitions is the presence of sharp interfaces at which the main thermophysical and optical characteristics: enthalpy H, heat capacity C_p and thermal conductivity λ, density ρ, pressure p, and reflectivity R of the surface undergo a steplike change.

In the case of rapid phase transformations that typical for pulsed nanosecond action, a complete set of equations for mass, momentum, and energy flows is used to describe heterogeneous phase transitions. However, considering that the absence of phase equilibrium at interphase boundary is a necessary condition for phase transition to occur, the expressions for conservation laws must be supplemented by appropriate kinetic relations characterizing the degree of nonequilibrium of phase transition.

2.1.1. Left fixed boundary, x = Г_s

The condition that flow of mass and heat be equal to zero is used as boundary conditions on the left fixed boundary:

$x={\Gamma }_s:u=0,W_T=0;$ E2

2.1.2. Model of heterogeneous (surface) melting: crystallization, õ = Г_sℓ (t)

For fast phase transitions, heterogeneous melting model consists of a system of equations expressing three conservation laws: mass, momentum, and energy supplemented by kinetic condition for the velocity of melting front υ_sℓ [57] obtained from molecular-kinetic theory [58] and which is the main characteristic of the process melting-crystallization. In a stationary (laboratory) coordinate system, surface melting-crystallization model written on moving interface of melting x = Г_sℓ (t) can be represented as:

$\begin{array}{l}x={\Gamma }_{\mathit{s\ell }}\left(t\right):{\rho }_s\left(u_s-{\upsilon }_{\mathit{s\ell }}\right)={\rho }_{\ell }\left(u_{\ell }-{\upsilon }_{\mathit{s\ell }}\right)\\ p_s+{\rho }_s{\left(u_s-{\upsilon }_{\mathit{s\ell }}\right)}^2=p_{\ell }+{\rho }_{\ell }{\left(u_{\ell }-{\upsilon }_{\mathit{s\ell }}\right)}^2\\ {\left(\lambda \left(T\right)\frac{\partial T}{\partial x}\right)}_s-{\left(\lambda \left(T\right)\frac{\partial T}{\partial x}\right)}_{\ell }={\rho }_s{L}_m^{\mathit{ne}}{\upsilon }_{\mathit{s\ell }}\end{array}$ E3

${\upsilon }_{\mathit{s\ell }}\left(\Delta T_{\mathit{s\ell }}\right)=\alpha {\left(3k_B{T}_{\mathit{s\ell }}/m\right)}^{1/2}\left(\text{exp}\left(\beta \frac{L_m}{k_B{T}_m\left({\mathrm{p}}_s\right)}\frac{\Delta T_{\mathit{s\ell }}}{T_{\mathit{s\ell }}}\right)-1\right)$ E4

$L_m^{\mathit{ne}}=L_m\left(T_m\left(p_s\right)\right)+\Delta C_{\mathit{ps}}\Delta T_{\mathit{s\ell }}+\frac{{\rho }_s+{\rho }_{\ell }}{{\rho }_s-{\rho }_{\ell }}\frac{{\left(u_s-u_{\ell }\right)}^2}{2},L_m\left(T_m\left(p_s\right)\right)=L_{m,0}+\delta \left(T_m\left(p_s\right)-T_{m,0}\right)$ , ΔC_ps  = C_ps  − C_pℓ  , ΔT_sℓ  = T_sℓ  − T_m (p_s ) , T_m (p_s ) = T _m , 0 + θp_s ,

where $L_m^{\mathit{ne}}$ is nonequilibrium heat of melting, L _m,0, T _m,0 are equilibrium heat of melting and melting point, respectively, p_s is pressure on a solid surface. α, β, and δ are parameters determined from molecular modeling [59, 60]. For Al, α = 0.21, β = 5.28, δ = 6.37 J mol^{− 1} K⁻¹, θ = 6.44 × 10⁻³ K/atm.

2.1.3. Model of heterogeneous evaporation x = Г_lv

Investigation of heterogeneous evaporation process began already in the century before last from experimental work: Hertz [61] and theoretical work: Knudsen [62] and continues at the present time, which is determined by its practical importance and not fully clarified features of nonequilibrium behavior of the substance during its evaporation. The mechanism of heterogeneous laser evaporation is realized in subcritical region of surface temperature T_sur and saturated vapor pressure p_sat = p_sat (T_sur ), for which the following inequalities hold: T_sur < T_cr and p_sat < p_cr . In the case T_sur > T_cr or p_sat > p_cr , then a supercritical laser evaporation regime is realized, at which the state of the substance varies continuously [63].

Heterogeneous evaporation is characterized by a sharp phase boundary. Thin (several mean path lengths) nonequilibrium Knudsen layer (KL) is adjacent to the boundary surface. The nonequilibrium of the KL is determined by the flow of matter through the phase boundary. The total flow consists of flow of particles (atoms, molecules) leaving the condensed phase and flow of particles (atoms of molecules) returning from the evaporated matter as a result of collisions. These flows have different distribution functions, f ⁽⁺⁾ and f ⁽⁻⁾ respectively, and that leads to a strong nonequilibrium. The flow of injected particles satisfies the Maxwellian distribution function f ⁽⁺⁾ with parameters of the surface of the condensed medium. The return flow is determined by conditions of gas-dynamic expansion outside the KL. The behavior of particles inside a nonequilibrium Knudsen layer is described by Boltzmann equation, by solving of which one can determine the unknown parameters T_υ and ρ_υ that are boundary conditions for the continuum equations of gaseous medium. In the heterogeneous evaporation model, a special role is played by parameter M (M = u_υ/u_sound is Mach number) on the outer side of KL, which determines the degree of nonequilibrium of the phase transition.

In the phase equilibrium state, when the pressure of saturated vapor p_sat (T_sur ) is equal to the external pressure p_υ , the parameter M = 0. In the subsonic evaporation regime, when M < 1, the behavior of the interface depends on gas-dynamic perturbations in the vaporized matter flow. The maximum nonequilibrium is determined by the maximum value of mass flow, which is known to be achieved at M = 1, when flow of matter through the interface is maximal, and the recoil pressure at the heated surface is minimal. In this case, such evaporation regime can be realized when the behavior of condensed medium no longer depend on the external gas-dynamic problem, which greatly simplifies the description of evaporation process. However, the implementation area for such regime remains uncertain, because due to the nonlinearity of the gas dynamics equations, the transition line between M < 1 and M ≥ 1 depends on desired solution and cannot be determined in advance [64]. Taking into account strong spatiotemporal diversity of processes in the nonequilibrium layer and continual media, KL is usually represented as a strong discontinuity in the gas-dynamic parameters. In such case, the kinetic processes in KL are not explicitly considered, and various phenomenological approaches [65–72] were used to determine boundary conditions on the outer side of the KL, which makes it possible without solving the kinetic problem, to determine the joining conditions under certain assumptions about the type of the nonequilibrium distribution function inside the discontinuity. One of the first works in which intensive evaporation was analyzed on the basis of phenomenological model is the work of Crout [65]. In this work, a nonequilibrium ellipsoidal Maxwellian particle distribution function was used written in an analytical form with the anisotropic in spatial directions longitudinal and perpendicular temperatures. Later papers used Mott-Smith approach [66], which was applied to the structure of the shock waves. With the help of this approach [67, 68] expressions for gas-dynamic parameters at M = 1 were obtained. Later in [69], this approach was extended to the entire range of evaporation 0 ≤ M ≤ 1. In [70, 71], the approximation of the distribution function over the entire evaporation range 0 ≤ M ≤ 1 was carried out by other more complicated expressions that satisfied the additional requirements that were formulated taking into account basic laws of gas dynamics. Following this laws, the used distribution functions must when M → 1 to provide extreme values for all three complete flows of mass j_m , momentum j_i and energy j_e . The model [65] also satisfies these requirements. At the same time, in widely used model [69], the requirement of extrema of all flows at selected point M = 1 is not fulfilled. The calculations showed that total flows, j_m , j_i , j_e depending on M have extrema at M = 0.88, 1.18, 1.22, respectively. Nonfulfillment of the requirement of extremes indicates an unsuccessful choice of distribution function for the reverse flow f ⁽⁻⁾. The mathematical model of surface evaporation in the Knudsen layer approximation consists of three conservation laws and two additional relations from which the parameters on outer side of the Knudsen layer are determined: temperature T_υ and density ρ_υ . Velocity u_υ and Mach number M are determined from the solution of the equations of gas dynamics:

$\begin{array}{l}x={\Gamma }_{\mathit{l\upsilon }}(t):j_{\mathit{\ell \upsilon }}^m={\rho }_{\ell }\left(u_{\ell }-{\upsilon }_{\mathit{\ell \upsilon }}\right)={\rho }_{\upsilon }\left(u_{\upsilon }-{\upsilon }_{\mathit{\ell \upsilon }}\right)\\ j_{\mathit{\ell \upsilon }}^i=p_{\ell }+j_{\mathit{\ell \upsilon }}^m\left(u_{\ell }-{\upsilon }_{\mathit{\ell \upsilon }}\right)=p_{\upsilon }+j_{\mathit{\ell \upsilon }}^m\left(u_{\upsilon }-{\upsilon }_{\mathit{\ell \upsilon }}\right)\\ j_{\mathit{\ell \upsilon }}^e=-j_{\ell }^T+j_{\mathit{\ell \upsilon }}^m\left[H_{\ell }+\frac{{\left(u_{\ell }-{\upsilon }_{\mathit{\ell \upsilon }}\right)}^2}{2}\right]-\sigma T^4=-j_{\upsilon }^T+j_{\mathit{\ell \upsilon }}^m\left[H_{\upsilon }+\frac{{\left(u_{\upsilon }-{\upsilon }_{\mathit{l\upsilon }}\right)}^2}{2}\right]\end{array}$ E5

where $j_{\ell }^T=W_{\ell }=-\lambda \left(T_{\ell }\right)\frac{\text{∂}{T}_{\ell }}{\text{∂}x},j_{\upsilon }^T=W_{\upsilon }=-\lambda \left(T_{\upsilon }\right)\frac{\text{∂}{T}_{\upsilon }}{\text{∂}x}$ . With these expressions taken into account, energy conservation law can be represented as:

$W_{\ell }-W_{\upsilon }={\rho }_{\ell }{\upsilon }_{\mathit{\ell \upsilon }}{L}_{\upsilon }^{\mathit{ne}}+\sigma T^4,$ E6

where $L_{\upsilon }^{\mathit{ne}}=L_{\upsilon }^e\left(T_{\ell }\right)+C_{\mathit{p\upsilon }}\left(T_b-T_{\mathit{\ell \upsilon }}\right)+\frac{{\rho }_{\ell }+{\rho }_{\upsilon }}{{\rho }_{\ell }-{\rho }_{\upsilon }}\frac{{\left(u_{\ell }-u_{\upsilon }\right)}^2}{2}$ is the nonequilibrium heat of evaporation, σ is Stefan-Boltzmann law constant.

To determine parameters T_υ , ρ_υ , p_υ , a modified Crout model was used [71].

$T_{\upsilon }={\alpha }_{{\rm T}}\left(M\right)T_{\mathit{sur}},{\rho }_{\upsilon }={\rho }_{\mathit{sat}}{\alpha }_{\rho }\left(M\right),$ E7

$\begin{array}{l}{p}_{\mathit{sat}}=p_b\text{exp}\left[\frac{L_{\upsilon }\left(T\right)}{k_B}\left(\frac{1}{T_b}-\frac{1}{T_{\mathit{sur}}}\right)\right],{\alpha }_{{\rm T}}\left(M\right)=\frac{2\gamma M^2{\left(m^2+0.5\right)}^2}{{\left(1+\gamma M^2\right)}^2{m}^2{t}^2},\\ {\alpha }_{\rho }\left(M\right)=\frac{1}{\text{exp}\left(-m^2\right)+{\pi }^{1/2}m\left(1+\text{erf}\left(m\right)\right)}\cdot \frac{\left(1+\gamma M^2\right)m^2}{\gamma M^2{\left(m^2+0.5\right)}^2},\end{array}$

the value m is determined from the equation F(M)(m ² + 0.5)² − m ²(m ² + 1.5 + a) = 0, where $F\left(M\right)=1+\frac{3\gamma M^2-1}{{\left(\gamma M^2-1\right)}^2},$ a = 2t²−0.5π^1/2mt−1, $t=\frac{2m}{{\pi }^{1/2}}+\frac{1+\text{erf}\left(m\right)}{\text{exp}\left(-m^2\right)+{\pi }^{1/2}m\left(1+\text{erf}\left(m\right)\right)}$ .

ρ_sat and p_sat are the saturated vapor density and pressure, p_b is the pressure under normal conditions (1.01325 × 10⁵ Pa), u_sound  = (γRT)^1/2 is the sound velocity, M = u_υ/u_sound is the Mach number, T_b is the boiling temperature, k_B is Boltzmann’s constant. Taking into account p_sat , the value of ρ_sat was determined from the equation p = p(ρ,T). For the laser radiation transfer equation on the target surface, the condition $G\left({\Gamma }_{\mathit{l\upsilon }}\right)=\left(1-R\left(T_{\mathit{sur}}\right)\right)G_0\text{exp}\left(-{\left(\frac{t}{\tau }\right)}^2\right)$ was used where R(T_sur ) is the temperature dependence of reflectivity of target surface.

2.1.4. Model of surface condensation

When the inequality p_sat < p is reached, the direction of phase transition changes, evaporation is replaced by condensation process on the target surface, to which the value M < 0 on outer side of the KL corresponds. Unlike evaporation, surface condensation can also occur in supersonic regime. According to [73] in subsonic condensation regime, gas dynamic quantities T_υ , M, and p_υ are related only by one two-parameter dependence p_υ/p_sat = F(T,M), where T = T_υ/T_sur and M are dimensionless parameters of temperature and velocity determined by the state of gas dynamic flow far from KL and are extrapolated from the gaseous medium to the discontinuity surface. The function F(T,M) is determined in advance from the solution of Boltzmann equation, the results of which are then tabulated [74]. The values of the function F(T,M) depend weakly on the parameter T and much stronger on the parameter M. Table values from [74] can be approximated by the expression [75]:

$p_{\upsilon }/p_{\mathit{sat}}=F\left(T,M\right)\approx F\left(M\right)=0.95\times \text{exp}\left(-2.42M\right)$ E8

When passing through the point M = −1 to supersonic regime of surface condensation, the boundary conditions change. In this case, all quantities on outer side of KL depend on the state of gas medium far from it and are extrapolated.

Thus, the description of the kinetics of heterogeneous evaporation is carried out by two one-parameter dependencies used as boundary conditions (for a known parameter M), and subsonic surface condensation is described by only one two-parameter (M,T) dependence.

2.1.5. Model of volumetric boiling of liquid phase heated by a laser pulse

The greatest difficulty in the continuum approach is the description of homogeneous mechanisms of phase transformations: melting crystallization and evaporation. Homogeneous mechanisms of phase transformations are characterized by nucleation of a new phase in a certain volume of superheated/supercooled matter. Representing them in continuum hydrodynamic models requires considerable additional efforts [57, 76] associated with formation of a cavity filled with vapor within a condensed medium.

The simplest scheme for simulating volumetric boiling in one-dimensional approximation can be represented by introducing, when certain criteria are satisfied into superheated liquid phase of artificial quasi-nuclei with thickness h_i (t) bounded by moving planes of liquid-vapor interface x_i (t) = Γ _{ℓυ , i} (t) ,  x _i + 1(t) = Γ _{ℓυ , i + 1}(t) where i = 1, 2, … is quasi-nuclei number. The criterion for nucleation beginning is a moment when superheating limit temperature T _ℓ,max of liquid phase is reached and spatial coordinate. First quasi-nucleus of vapor phase with initial width Δx _υ 1(t) = Γ _{ℓυ , 2}(t) − Γ _{ℓυ , 1}(t) ~ 5 nm is placed at this point. For a new region of vapor bounded by two flat surfaces, the initial conditions T _υ , 1 = T _{ℓ , max} ,  ρ _υ , i  = ρ_sat (T _{ℓ , max}) ,  p _υ , 1 = p_sat (T _{ℓ , max}) were set. Under the influence of pressure difference between the inside quasi-nuclei and on the external irradiated surface of the liquid layer with the temperature T_sur , a rapidly expanding cavity filled with vapor is formed and blowout the liquid layer with thickness Δx _ℓ 1(t) = Γ _{ℓυ , 1}(t) − Γ_ℓυ (t) in the direction of gaseous medium. The description of the processes in the expanding cavity and in the near-surface liquid layer separated from main target was carried out using gas-hydrodynamic system of Eq. (1). The model of heterogeneous evaporation (5)–(7) was used as boundary conditions on interphase planes Γ _{ℓυ , i} (t) ,  Γ _{ℓυ , i + 1}(t). Also, absorption and reflection of the intensity of laser radiation in the formed layer were taken into account. In case of repeated explosive boiling with the formation of next liquid fragment for each of them, as well as the main target, gaseous medium and formed vapor cavities, the solution algorithm remains unified.

The values of limit of superheat temperature depend on the rate of energy input and can be determined in advance from the molecular dynamics simulation. For Al, depending on the rate of heating, limit of superheat temperature of liquid phase T_max is in the range T _max ~ (0.89 ÷ 0.95)T_cr [77].

It should be noted that the use of gas dynamic equations inside the cavities which dimensions are small in the initial moments of time in comparison with the mean free path is a simplifying approximation that makes it easy to take into account the influence of expansion velocity of the cavity.

2.1.6. Moving contact boundary, x = Г_vg (t)

On contact boundary and front of the shock wave, well-known standard relations are used [78]. At vapor-air interface, the boundary conditions were set in the form of equal values of velocity, pressure, and temperature:

$x={\Gamma }_{l\upsilon }\left(t\right):u_{\upsilon }=u_g,P_{\upsilon }=P_g,T_{\upsilon }=T_g$

2.1.7. Moving shock wave, x = Г _sh,g (t)

A shock wave in air Г_sh,g (t) is a strong nonstationary discontinuity, on which three conservation laws are written in laboratory coordinate system [78]:

$\begin{array}{l}{j}_{\mathit{sh},g}^m={\rho }_1\left(u_1-{\upsilon }_{\mathit{sh},g}\right)={\rho }_0\left(u_0-{\upsilon }_{\mathit{sh},g}\right),\\ j_{\mathit{sh},g}^i=p_1+{\rho }_1{\left(u_1-{\upsilon }_{\mathit{sh},g}\right)}^2=P_0+{\rho }_0{\left(u_0-{\upsilon }_{\mathit{sh},g}\right)}^2,\\ j_{\mathit{sh},g}^m\left({\epsilon }_1+\frac{{\left(u_1-{\upsilon }_{\mathit{sh},g}\right)}^2}{2}\right)-W_{T,1}=j_{\mathit{sh},g}^m\left({\epsilon }_0+\frac{{\left(u_0-{\upsilon }_{\mathit{sh},g}\right)}^2}{2}\right)-W_{T,0}\end{array}$ E9

The indices 0 and 1 denote the values of the quantities on side of the background and the shock wave, respectively.

2.1.8. Right moving boundary, x = Г_g (t)

The right boundary on the side of unperturbed gas is declared moving in order to improve economic effectiveness of the computational algorithm [79]. The speed of its motion is found from the differential equation of momentum.

2.2. Computational algorithm

The differential model (1)–(8) was approximated by a family of conservative finite-difference schemes [80, 81] written on computational grids with dynamic adaptation [79, 82, 83]. The method is based on the idea of transition to an arbitrary nonstationary coordinate system that allows calculations with an arbitrary number of discontinuous solutions, such as shock waves, propagating phase and temperature fronts, contact boundaries, and spalled fragments.

3. Temperature dependences of aluminum properties

In the mathematical modeling of the process of laser action on metals with an absorption coefficient of κ_L  ≥ 10⁷ m⁻¹, the model of surface heating and evaporation is widely used. In this model, the maximum of the temperature profile T_max coincides with surface temperature T_sur from which evaporation occurs. At lower values of κ_L typical for dielectrics as well as for metals when irradiated with an electron beam or X-ray pulsed radiation or in the case of a metal-dielectric transition [84, 85], bulk heating of the substance and surface evaporation leads to the formation of temperature maximum below the surface of the substance. As a result, a volumetric explosive boiling of superheated liquid can occur in the region of the temperature maximum. The explosive boiling of superheated liquid is closely related to the concept of metal-dielectric phase transition. Zeldovich and Landau [86] denoted the possibility of a metal-dielectric phase transition for an expanded metal at subcritical temperatures. To perform the metal-dielectric transition, it is necessary that the energy acquired by metal atoms is higher than their binding energy in the crystal lattice and the distance between atoms is equal to the value causing a violation of their short-range order, which leads to the localization of the electrons on atoms [87].

With the advent of lasers capable of evaporating metals and producing plasma on their surfaces, it has become possible to observe laser-induced phase transitions metal-dielectric with the formation of transparency waves in the massive targets [88] and thin metallic films [89] in the subcritical temperature range. In other experimental and theoretical studies of the interaction of ns-laser pulses with an intensity of 10⁷–10⁸ W/cm² with metallic targets, the results of the appearance of metal-dielectric transitions accompanied by the formation of transparency waves are reported [53, 54, 85, 88–91].

Nevertheless, the physical mechanisms of interaction of laser radiation with metals taking into account metal-insulator transition have not yet been fully studied. This causes great difficulties in determining thermophysical and optical characteristics of metals in the vicinity of the critical region.

Figure 2 shows temperature dependences of thermophysical λ(T), L_v (T), C_p (T) and optical R(T), κ_L (T) characteristics of aluminum. The curves were L_v (T) and C_p (T) obtained from molecular dynamics calculations, and λ(T), κ_L (T), and R(T) were constructed on the basis of theoretical concepts [92–94] and reference data [95–97].

Figure 2.

Temperature dependences of Al properties (a–d).

4. Results discussion

One of the purposes of this work is a detailed study of mechanism of explosive boiling in metals since explosive boiling is considered to be the most efficient thermal mechanism for the laser ablation of materials. At the same time, the difficulties associated with understanding of the mechanism of homogeneous phase transitions in metals occurring under the action of ns-laser pulses are known. In [98], based on analytical solution of thermal model, it was established that laser removal of material from a solid target, with a certain choice of irradiation parameters and material, is determined by the presence of near-surface temperature maximum. The validity of this position was confirmed by the results of numerical solution of thermal model for low-absorbing liquids irradiated by laser pulses [43] and nonmetallic solid materials [47]. For strongly absorbing media mainly metallic, calculations based on thermal model [44, 45, 47, 56] have shown that the magnitude of near-surface temperature maximum is several degrees. On this basis, overheating was excluded from consideration, up to the statement [56] that in metals, maximum temperature is always on the surface of the target, and sub-surface superheating is impossible. Explosive boiling was interpreted as a surface spitting of liquid phase when a critical temperature is reached on the surface [99]. This interpretation is not convincing, since it does not allow to determine even approximately the parameters of explosive boiling.

4.1. Mode of exposure

Let us consider ns-PLA process of aluminum target under the mode of exposure close to experimental conditions [53]. Al target with the thickness of 10⁻⁴ m is placed in the air under the normal conditions with pressure 1 bar and room temperature (T = 300 K). Laser pulse of Gaussian shape G = (β/π)^1/2 G ₀ exp(−β (t/τ)²) with full width at half maximum (FWHM) τ = 5 × 10⁻⁹ s, with wavelength λ_L  = 1.06 μm and fluence F = 3.5 J/cm², where G ₀ = 6.1 × 10⁸ W/cm² is the peak intensity at t = 0, − ∞  < t < ∞, β = 4 × ln2 falls on the target surface from right to left. Air for the selected exposure mode is completely transparent for laser radiation, which is partially reflected from metal surface and partially absorbed by target material layer. The energy release of laser pulse has volumetric nature.

4.2. Evolution of processes

At the initial stage, time evolution of the processes at the surface of the target, in the target and in the gas medium near the target such as appearance of phase fronts (melting, evaporation), contact boundary, and shock wave and associated with formation of new phase media (liquid, vapor) occurs at the leading edge of laser pulse. Heterogeneous melting begins at the moment t = −3.9 × 10⁻⁹ s. The maximum velocity of melting front reaches a value υ_sl  = 130 ms⁻¹. As the interior of the target is heated, the melting front Г_sl (t) runs from its surface forming a new region of liquid phase. Further heating leads to appearance of heterogeneous evaporation front Г_lv (t) that runs inside the melt. A flow of the evaporated matter was formed at the surface of the melt pushing out air and creating another new phase—vapor. The new area occupied by vapor is limited on the one hand by the moving interface Г_lv (t) (evaporating surface) and on the other hand by the moving contact boundary vapor-air Г_vg (t). The surface evaporation process is controlled by the surface temperature T_sur (t) and Mach number M(t) on the outer side of KL (0 < M(t) ≤ 1). Evaporation begins at the moment t = −2.3 × 10⁻⁹ s when the saturated vapor pressure exceeds the pressure of external gas p_sat (t) > p_g (t). The maximum velocity of evaporation front υ_ℓv  = 88 ms⁻¹ is approximately 1.5 times smaller than velocity of melting front. The vaporized matter flow acting as a piston pushes out cold air, and performing a certain work warms up at peak G ₀ intensity up to temperature T_vap  = 4.7 × 10³ K. Under the pushing action of vapor flow, the compression of cold dense air which turns into a shock wave occurs, t = − 2.1 × 10⁻⁹ s. The spatial structure of erosion plume is shown in Figure 3 by the main characteristics T(x), ρ(x), u(x), p(x) after the formation of shock wave.

Figure 3.

The spatial profiles of (a) T(x), (b) ρ(x), (c) u(x), and (d) p(x) in the target and the gas medium at the moment t = −1.9 ns.

By the moment t = 0, the shock wave propagates with velocity υ_sh,g  = 2.9 km s⁻¹ and temperature T_sh,g  = 4.4 × 10³ K, ahead of contact boundary moving with speed υ_sh,g  = 2.5 km s⁻¹ toward the laser flow. The temperatures of vapor and air are in this case insufficient for initiation of ionization, and vapor-gas medium remains transparent for laser radiation.

4.3. Formation of near-surface temperature maximum

A consistent study of laser ablation is complicated by the fact that heterogeneous and homogeneous mechanisms of both melting and evaporation prove to be interrelated, and this interaction must be taken into account explicitly. Pulsed laser action on materials (including metals) has volumetric nature of energy release. Thus, when target is heated strongly in the vicinity of the critical point (T_cr  = 7600 K, ρ_cr  = 0.47 g cm⁻³, p_cr  = 1.42 kbar), thermophysical and optical properties of liquid phase change abruptly. Heat capacity C_p when approaching a critical point tends to infinity C_p → ∞, thermal conductivity λl and absorption κ_L coefficients tend to zero, surface absorptivity A = (1−R) → 1, Figure 2a – c . Under the joint influence of volumetric energy release of laser pulse and energy transfer by surface flow of evaporating substance, a near-surface maximum of temperature T_ℓ,max is formed in the depth of liquid phase, for which the following relations are satisfied: T_ℓ,max > T_sur , p_ℓ,max,sat > p_sur,sat.

Figure 4a shows a fragment of the spatial temperature distribution in near-surface layer of the target with overheating T_ℓ,max − T_sur  = 170 K at a depth of 30 nm from the irradiated surface at the moment immediately preceding the explosive boiling. For comparison, Figure 4b shows a fragment of temperature profile, calculated under the same conditions with help of molecular dynamics, which showed close results. The maximum was located at a depth of 70 nm, and the superheat value was 150 K. Thus, simulation results obtained by different methods indicate the presence of temperature inhomogeneity in metallic target (Al) caused by a decrease in the temperature of the irradiated surface by heterogeneous evaporation. Thus, the conditions for phase explosion were created with formation of a cavity in the region of maximum temperature, where the deepest passing into the region of superheated metastable liquid critical point of liquid-vapor transition is achieved.

Figure 4.

Spatial profiles of temperature maximum, (fragments) before explosive boiling: (a) continual model and (b) atomistic model.

4.4. Explosive boiling

The process of explosive boiling begins at the backside front of the laser pulse. First boiling occurs at the moment t = +1.5 ns when the maximum of the temperature profile reaches the temperature of limiting superheating of liquid phase equal to T_ℓ,max  = 0.9T_cr  = 6840 K. Figure 5a – d shows the spatial profiles of T(x), ρ(x), u(x), p(x). Arising cavity is the result of homogeneous nucleation in superheated liquid phase and is accompanied by sharp increase in pressure p_max = p_sat (T_ℓ,max ) which is more than two times greater than the recoil pressure on irradiated free surface p_sur  = 0.55p_sat (T_sur ). This pressure plays a decisive role in the expansion of cavity filled with vapor and in rapid growth of flyout velocity of thin d ₁ = 30 nm fragment of liquid metal. Because of small thickness, flying fragment weakly absorbs laser radiation, the main part of which is released in the target increasing its temperature. Both surfaces of spalled fragment and surface of the target are subject to intense surface evaporation. Heterogeneous evaporation leads to decrease in temperature of the fragment and causes strong cooling of target surface, contributing to formation of next near-surface temperature maximum in which a new explosive boiling occurs when the temperature of limit overheating T_ℓ,max  = 0.9T_cr of metastable liquid is reached. Repeated boiling occurred at a time t = +1.7 ns. The total number of explosive boilings reaches 5. The moments of occurrence and thickness of outgoing fragments are given in Table 1 .

Figure 5.

The spatial profiles of (a) T(x), (b) ρ(x), (c) u(x), (d) p(x) in the target and the gas medium at the moment t = −1.7 ns.

Fragment number, n	1	2	3	4	5
Moment of occurrence, ns	1.5	1.7	1.9	2.1	2.4
Fragment thickness d, nm	30	37	26	21	20

Table 1.

Number, time of appearance and thickness of spalled fragments.

Figure 6a – d shows the state of all liquid fragments of explosive boiling at the time t = +3.0 ns. The resulting fragments eventually acquire a sufficiently high flyout velocity and under the influence of laser radiation absorption continue to evaporate until they disappear completely. The general trend in the evolution of all fragments is that over time the cavity size, flyout velocity and density of matter increase, and pressure in the cavity, thickness and temperature of the fragments decrease. So in the initially formed cavity the width increased from 5.0 nm to 1.7 μm, the velocity increased up to 2.3 km/s and density due to cooling changed from 0.65 g/cm³ to 1.3 g/cm³. At the same time, the pressure in the cavity decreased to 200 bar, the thickness of fragment on the verge of extinction due to surface evaporation fell to h_ℓ  = 1.9 nm, and the temperature decreased to T_ℓ  = 5400 K.

Figure 6.

The spatial profiles of (a) T(x), (b) ρ(x), (c) u(x), (d) p(x) in the target and the gas medium after five explosive boiling at the moment t = −3.0 ns.

The simulation results are in good qualitative agreement with the experimental data [53]. Quantitatively, the total depth due to release of superheated liquid phase (d_Σ  = 37.0 nm) differs by about two times in excess. The total removal of matter by the mechanism of explosive boiling considerably exceeds the amount of evaporated matter. The trajectory of shock wave motion almost completely coincides with the experimental dependence [53] up to the moment t = 200 ns. From now, calculated curve is located above the experimental curve. As the fluence increases up to F ~ 10 J/cm², the target is heated faster and explosive boiling process begins to move to leading edge of the laser pulse. The frequency and amount of boiling thus increase to tens and hundreds of ejected liquid fragments, but the overall picture of processes does not undergo qualitative changes.

5. Conclusion

Within the framework of continuum approach, a multiphase, multifront hydrodynamic model and computational algorithm were constructed that were designed for modeling ns-PLA of metal targets embedded in gaseous media. The model contains temperature dependences of optical, thermal, and thermodynamic properties of metal target (Al). The temperature dependences of heat capacity C_p (T) and thermal conductivity λ(T) of aluminum take into account their singularity in the vicinity of critical point. In optical characteristics, the metal-insulator phase transition was taken into account accompanied by the formation of the transparency wave.

The computational algorithm is based on the method of dynamic adaptation, which ideally fit to the study of problems with heterogeneous phase transitions and allows the production of numerical solutions with an explicit allocation of unlimited number of interface boundaries in condensed medium, contact boundaries, and shock waves in a gaseous medium.

The calculations have shown that under the influence of temperature dependences of optical and thermophysical properties of Al the nature of laser heating of the target in high-temperature region changes substantially. The energy release has a volume nature, and irradiated surface is markedly cooled by the process of surface evaporation. An inhomogeneity arises in the spatial temperature profile—a near-surface maximum of temperature. When the limiting superheating temperature is reached, the conditions for explosive boiling were realized with formation of a cavity in the region of the maximum, where the maximum penetration into the region of superheated metastable liquid near the critical point of liquid-vapor transition is achieved. For Al in action modes under the consideration, the typical depths of temperature maximum are d = 20–70 nm with the values ΔT_max = T_ℓ,max − T_sur = 150–170 K.

An approach was proposed to formulation in continuum models of homogeneous nucleation in a superheated liquid phase based on the introduction into hydrodynamic model of quasi-nuclei with thickness of 1–5 nm. Generation of quasi-nuclei is carried out by the criterion of maximum permissible overheating of initial region. The criterion of limit overheating is determined from molecular dynamics simulation.

Mathematical modeling using the developed technique allowed to obtain a sequence of five explosive boilings for the threshold fluence F = 5.2 J/cm² and 14 for F = 6 J/cm². In the first case, the depth of ablation due to explosive boiling was 37 nm, in the second—70 nm. Due to the surface evaporation, they were 17 nm and 12 nm, respectively. The obtained data make it possible to conclude that phase explosion is the main mechanism of material removal in nanosecond range of laser action with fluence F ≥ 5.2 J/cm². The upper values of τ and F are limited by the processes of supercritical expansion and plasma formation and are subject to determination.

In metals, as in nonmetals, the transition from surface evaporation to volume removal of mass occurs in the region of temperature maximum near the critical temperature, where the influence of overheated metastable liquid phase determines the competition between surface evaporation and explosive boiling.

It should be noted that in explosive boiling regimes under investigation, cavity formation takes place in the region of positive values of pressure generated by homogeneous nucleation in superheated liquid phase, in contrast to the case of a series of spalles [57] arising in the region of negative pressures in unloading wave upon action to metal targets of ultrashort pico-femtosecond laser pulses.

Validation of obtained results was performed by comparison with experimental data and with calculations with atomistic models [77, 100].

Acknowledgments

The research was funded by the Russian Scientific Foundation (grant No 15-11-00032).