\documentclass[reqno]{amsart}
%\usepackage[notref,notcite]{showkeys}
\usepackage{hyperref}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small
Eighth Mississippi State - UAB Conference on Differential Equations and
Computational Simulations.
{\em Electronic Journal of Differential Equations},
Conference 19 (2010), pp. 1--14.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2010 Texas State University - San Marcos.}
\vspace{9mm}}
\begin{document} \setcounter{page}{1}
\title[\hfilneg EJDE-2010/Conf/19/\hfil Enthalpy model]
{Enthalpy model for heating, melting, and vaporization in laser ablation}
\author[V. Alexiades, D. Autrique \hfil EJDE/Conf/19 \hfilneg]
{Vasilios Alexiades, David Autrique}
\address{Vasilios Alexiades \newline
Mathematics, University of Tennessee,
Knoxville, TN 37996, USA \newline
\emph{and} Oak Ridge National Laboratory, Oak Ridge TN 37831, USA}
\email{alexiades@utk.edu}
\address{David Autrique \newline
Department of Chemistry, University of Antwerp,
Antwerp, Belgium}
\email{david.autrique@ua.ac.be}
\thanks{Published September 25, 2010.}
\subjclass[2000]{92C45, 35K60, 65M99}
\keywords{Laser ablation; phase change; enthalpy scheme; thermochemistry;
\hfill\break\indent
finite volume scheme}
\begin{abstract}
Laser ablation is used in a growing number of applications in various
areas including medicine, archaeology, chemistry, environmental and
materials sciences.
In this work the heat transfer and phase change phenomena during
nanosecond laser ablation of a copper (Cu) target in a helium (He)
background gas at atmospheric pressure are presented.
An enthalpy model is outlined, which accounts for heating, melting, and
vaporization of the target.
As far as we know, this is the first model that connects the
thermodynamics and underlying kinetics of this challenging phase change
problem in a self-consistent way.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\section{Introduction}\label{S:1}
Laser ablation (LA) is the process of removing material from a target
using a laser beam. It is used in a growing number of applications,
such as pulsed laser deposition, nanoparticle manufacturing,
micromachining, surgery, as well as chemical analysis
\cite{ready71,vonAllmen87,chrisey94}.
Direct solid microanalysis using LA in combination with inductively
coupled plasma mass spectrometry (LA-ICP-MS) is nowadays one of the
most frequently used techniques for a fast and powerful multi-element
determination of solid samples at the trace and ultra-trace
concentration levels of a wide variety of sample types
\cite{vertes93,broekaert02,becker02,gunther02}.
The growing interest in LA as a sample introduction
technique stems from the ability to sample diverse materials,
ranging from conducting to non-conducting, inorganic and organic
compounds as solids or powders. Besides bulk analysis, the focussing
characteristics of lasers permit sampling in small areas, so that
localized microanalysis and spatially resolved studies
become feasible.
The processes occuring during LA-ICP-MS are depicted schematically in
Fig. \ref{fig:1}.
\begin{figure}[ht]
\centering
\includegraphics[width=0.95\textwidth]{fig1} % height=2.3in
\label{fig:1} \caption{
Schematic of the Laser Ablation - Inductively Coupled Plasma
Mass Spectrometry (LA-ICP-MS) process. (Courtesy of C.C. Garcia,
ISAS, Dortmund)}
\end{figure}
Within an ablation cell, an amount of matter is extracted from the
sample by focusing a laser beam on its surface.
The sample heats up, melts, vaporizes, and expands against
a carrier gas at atmospheric pressure.
The expanding plume partially ionizes under the laser radiation and a
plasma is generated above the target.
Complex physical processes occuring at the target surface and in the
expanding plume result in the formation of small particles with sizes
ranging from several nanometers to a few micrometers.
The aerosol is subsequently transported by the carrier gas through a
tubing system to an ICP-torch where it is ionized.
Finally a small fraction of the ions is sent to the MS-instrumentation
where they are analyzed.
The laser ablation process is typically investigated by means of
kinetic
\cite{noorbatcha87,ellegaard99,itina97, garrelie98,kools93},
hydrodynamic
\cite{bogaerts03,chen05,gusarov00,bleiner05,lindner09},
or hybrid models \cite{itina02,itina03}.
Because kinetic and hybrid models require quite long calculation times,
one tries to apply hydrodynamic models if at all possible.
Since the laser-target heating during LA-ICP-MS results in a cascade of
complex physical processes
(see \S\ref{S:2}),
it is of utmost importance to include in the existing hydrodynamical
models a consistent treatment of target heating.
Here the kinetic relations that interconnect the physical processes in
and above the target determine the pathway the metal follows on the
phase diagram during laser heating. Accordingly, one has to account
for them in order to arrive at a physically consistent treatment of
laser-material interaction.
In \S\ref{S:2} we briefly describe some issues that arise during the
modeling of \emph{nanosecond laser ablation}.
In \S\ref{S:3} we describe a seamless enthalpy-based model treating
the heating, melting and vaporization of the target.
Numerical simulations for laser ablation of copper (Cu) in a background gas (He) in 1-D, are presented in \S\ref{S:4},
and final conclusions in \S\ref{S:5}.
\section{Nanosecond laser ablation}\label{S:2}
Several complex, tightly coupled physical processes occur in and just
above the target. The target heats up, melts and vaporizes.
Above the target, the evaporated particles achieve translational
equilibrium within a few mean free paths by means of collisions.
This transition layer is known as the \emph{Knudsen layer} (KL).
The vapor near the target is a partially ionized gas that is separated
from the surface by an electrical sheath across which there is a
potential drop that accelerates the charged particles.
The Knudsen layer provides the connection between the `inner world'
(target) and the `outer world' (plasma plume).
Beyond this layer, the plasma absorbs laser energy before it can reach
the target, giving rise to plasma shielding.
The absorption of laser energy results in very high plume temperatures,
velocities, species densities, pressures, and so on.
At later times, small particulates are formed during the cooling stage
of the plume. Here homogenous nucleation and recondensation result in
nanosized particles. Moreover, recoil pressures acting on the melt can
cause melt motion and melt expulsion forming larger particles.
The melt pool dynamics leads to formation of jets at the edge of the
formed crater, that break up in micro-droplets
\cite{vonAllmen87,bleiner06}.
An effective hydrodynamic model of the entire ablation process
should posess the following highly desirable features:
\begin{itemize}
\item seamless treatment of heating / melting / vaporization,
\item temperature-dependent thermophysical and optical properties,
\item treatment of vapor as a non-ideal gas by means of an equation of
state(EOS) up to and above the critical temperature $T_{\rm crit}~$,
\item treatment of the plasma-sheet at the target surface,
\item coupling target and plume by means of the Knudsen layer,
\item treatment of plume expansion and recondensation, in vacuum or in a background gas
\item dealing with plasma formation (ions and electrons) mechanisms,
\item treating absorption of laser energy in the plume and shielding of the target,
\item capturing (very) strong shocks via high resolution numerical schemes,
\item be applicable in 1D, 2D axisymmetric and 3D spatial dimensions.
\end{itemize}
\noindent
Additional challenges in modeling laser ablation arise from:
\begin{itemize}
\item extreme space and time scales,
\item extreme gradients: temperature may rise to thousands of degrees locally,
\item extreme variation in thermophysical properties,
\item the need for extensive thermophysical data:
T-dependent density, heat capacity, thermal conductivity,
phase diagram for solid, liquid, vapor
over the range 300 K to $T_{\rm crit}$ ( $\sim 8000$ K ),
\item the need for T-dependent and wavelength dependent optical data.
\end{itemize}
\section{Enthalpy model for target heating - melting - vaporization}\label{S:3}
In this section we develop an enthalpy formulation for
target heating - melting - vaporization which addresses the issues
listed above.
Heat conduction and phase changes in the target are naturally
described by the energy conservation law in terms of the (volumetric)
enthalpy $H$.
\begin{equation}
\partial{}_t H + \nabla \cdot \vec{Q} = S \, .
\end{equation}\label{eq:3:1}
%
The laser source term $S$, at time $t$ and position $z$, is given by
%
\begin{equation}
S(t, z) = I(t) (1-R) \, \alpha e^{- \alpha z} ~,
\end{equation}\label{eq:3:2}
%
with $R$ and $\alpha$ the reflection and absorption coefficients of the
target, respectively, and intensity
%
\begin{equation}
I(t) = \beta I_o \, e^{ - (t-t_{max})^2 / 2\sigma^2 }\, .
\end{equation}\label{eq:3:3}
%
Here $\beta$ is the shielding coefficient at the target that depends
on the plume temperature and species densities
(neutrals, ions and electrons) above the target.
The calculations presented later in \S\ref{S:4} are performed for a
Gaussian-shaped laser pulse of full width at half maximum (FWHM) 10 ns,
and peak laser intensity of the order of $10^{13} {\mbox W/m^2}$,
at a wavelength of 266 nm.
In Fig. \ref{fig:2} the laser intensities before ($\beta = 1$) and
after shielding at the surface are shown.
%
\begin{figure}[ht]
\centering
\includegraphics[width=0.95\textwidth]{fig2} %,height=4.0in
\label{fig:2}\caption{
Time course $I(t)$ of laser pulse intensity with FWHM=10 ns,
peak $7 \times 10^{12} \mbox W/m^2$.
}
\end{figure}
There are three time scales that characterize laser heating:
$\tau_e = C_e / \gamma$, the electron cooling time,
$\tau_i=C_i / \gamma$, the lattice heating time, and
$\tau_l$, the duration of laser pulse.
Here $C_e$ and $C_i$ are the heat capacities of the lattice and
the electron subsystems, respectively, and
$\gamma$ a parameter characterizing the electron-lattice coupling in
the target. Since $\tau_e$ is of the order of 1 ps and
$\tau_l$ is of the order of ns, the electron absorbed energy has
enough time to be transferred to the lattice
\cite{chichkov96}.
Consequently, the electrons and the lattice reach thermal equilibrium
and heat conduction in the target can be described by means of
Fourier's law provided that the heat flux is constant over several
mean free paths of the electrons that carry the heat
\cite{harrington67}.
Since the melt thickness is very small, convective effects do not arise
in the melt. Thus the heat flux in the solid and liquid is due to heat
conduction only, $\vec{Q} = - k \nabla T~$,
with $k$ = thermal conductivity, $T$ = temperature.
On the contrary, the presence of strong advection in the plume domain
requires solution of the continuity and momentum equations as well.
Note that, depending on the local pressure $P$, vaporization can occur
at any temperature above absolute zero,
even below the melt point $T_m$ as well as above the normal
boiling point $T_b$.
Therefore the vaporization temperature $T_V (P)$ is unknown and must be
determined in the course of the solution.
Finally note that the surface recession velocity, which determines
the crater depth in the target, is also an unknown, since the surface
is an unknown moving boundary.
The energy conservation law \eqref{eq:3:1} is valid mathematically,
in the weak sense, throughout the system irrespective of phase
(with coefficients appropriate to phase:
$k^S$, $k^L$, $k^V$, etc.).
It directly updates the enthalpy $H$ (per unit volume) to new times,
from which we need to find the new phase, temperature, and
phase fractions. This requires Equations of State (EOS)
$H = H(T, P, phase)$ consistent with the thermochemistry of the
material, which we describe below.
\subsection{Phases from volumetric enthalpy}\label{S:3:1}
The development is based on the basic Gibbs relation for the exact
differential of enthalpy \cite{book} in any phase $i = S, L, V$ ~:
%
\begin{equation}
d H^i = {C_p}^i dT + [ -{\alpha}^i T + 1 ] dP \quad i=S, L, V~,
\end{equation}\label{eq:3:4}
%
where
$C_p^i$ and $\alpha^i$ are the heat capacity (per unit volume) and
thermal expansion coefficient of phase $i$, respectively.
The enthalpy within each phase (Solid, Liquid, Vapor) can be expressed
in terms of the pair ($T, P$) by integration along constant $T$ and
constant $P$ paths on the ($T, P$)-phase diagram,
of the Gibbs relation \eqref{eq:3:4}, relative to a \emph{consistent}
reference state.
In what follows, upper case $H$ denotes volumetric enthalpy
(per unit volume), lower case $h$ denotes specific enthalpy
(per unit mass), and values with subscript "${}_{\mbox o}$"
denote interfacial (switch) values.
Clearly, $H = \rho h$ in each phase.
We are dealing with three phases: solid (S), liquid (L), vapor (V).
The amount of phase $i = S, L, V$ present can be quantified by its
\emph{volume} phase fraction, $\Lambda^i$,
or by its \emph{mass} phase fraction, $\lambda^i$.
The fractions are related by: $\rho^i \Lambda^i = \rho \lambda^i$
($i = S,L,V$), where $\rho$ = total density (of mixture) =
$\Sigma_i \rho^i \Lambda^i$. Note that fractions must add up to 1:
$\Lambda^S + \Lambda^L + \Lambda^V = 1$ and
$\lambda^S + \lambda^L + \lambda^V = 1$.
%
For simplicity, we assume the density of solid, $\rho^S$, and of
liquid, $\rho^L$, remain constant (but $\rho^S \ne \rho^L$).
As reference state we choose \emph{solid} at temperature $T_m$ and
pressure $P_{\rm ref}=1$ atm. The value $h_{\rm ref}$
(and thus also $H_{\rm ref} = \rho^S h_{\rm ref}$) of the enthalpy
at the reference state can be arbitrary; it will be chosen later
(so as to make $\Delta H^{vap} > 0$).
By integration of \eqref{eq:3:4} within each phase, and changing phase
at ($T_m , P_{\rm ref}$) (from solid to liquid)
and at ($T_V , P_{\rm ref}$) (from liquid to vapor)
we can characterize the phases in terms of the value of $H$ alone
(which is being updated by the energy conservation law).
This is done in terms of certain {\bf\emph{switch}} values
$H_o^S$, $H_o^L$, $H_o^{LV}(T_V)$, and $H_o^V (T_V)$,
defined below, with
$T_V$ any desired vaporization temperature ($T_m < T_V \leq T_b$).
\vspace{2mm}
\noindent
$\bullet$ \textbf{Solid phase}\quad ($T \leq T_m$): \quad
$H < {\bf H_o^S} := H_{\rm ref} \equiv \rho^S h_o^S$, (set $h_o^S := h_{\rm ref}$).
\begin{center}
Then $\Lambda^S = 1$, $\Lambda^L = 0$, $\Lambda^V = 0$,
$\rho = \rho^S$, \quad
$H(T) = H_o^S + \int_{T_m}^T {C_p^{S} (\tau)} d\tau $. \\
Thus, given $H$, we can find the temperature $T$ ($< T_m$) by solving
the equation:\\
$\int_{T_m}^T {C_p^{S} (\tau)} d\tau = H - H_o^S$ \quad for $T$.
\end{center}
\noindent
$\bullet$ \textbf{Solid-Liquid transition}\quad (at $T = T_m$):\quad
${\bf H_o^S} < H < {\bf H_o^L} := \rho^L h_o^L$, \\
with $h_o^L = h_o^S + \Delta h^{fus}$. Set
$\Delta H^{fus} := H_o^L - H_o^S \equiv (\rho^L-\rho^S ) h_o^S + \rho^L \Delta h^{fus}$,
which can be made $> 0$ by choosing $h_{\rm ref}$ to be $\leq 0$.
%
\begin{center}
Then $T = T_m$, \quad $\Lambda^L = (H-H_o^L)/\Delta H^{fus}$,
$\Lambda^S = 1-\Lambda^L$, $\Lambda^V = 0$,
$\rho = \Lambda^S \rho^S + \Lambda^L \rho^L$,
$H = H_o^L + \Lambda^L \Delta H^{fus}$, $h = H/\rho$.
\end{center}
\noindent
$\bullet$ \textbf{Liquid phase}\quad ($T_m < T < T_V$):
${\bf H_o^L} < H < {\bf H_o^{LV}(T_V)} := H_o^L
+ \int_{T_m}^{T_V} {C_p^{L} (\tau)} d\tau $.
%
\begin{center}
Then $\Lambda^L = 1$, $\Lambda^S = 0$, $\Lambda^V = 0$,
$\rho = \rho^L$,
$H(T) = H_o^L + \int_{T_m}^T {C_p^{L} (\tau)} d\tau $. \\
Thus, given $H$, we can find the temperature $T$ ($T_m < T < T_V$)
by solving for $T$ the equation:\quad
$\int_{T_m}^T {C_p^{L} (\tau)} d\tau = H - H_o^L$.
\end{center}
\noindent
$\bullet$ \textbf{Liquid-Vapor transition}\quad (at $T = T_V$):\quad
${\bf H_o^{LV}(T_V)} < H < {\bf H_o^V (T_V, P_{\rm ref}) } :=\rho^V(T_V, P_{\rm ref}) h_o^V$,
with $h_o^V (T_V, P_{\rm ref}) := h_o^{LV}(T_V)
+ \Delta h^{vap}(T_V, P_{\rm ref})$, where \\
$\rho^V(T_V, P_{\rm ref})$ = vapor density at ($T_V, P_{\rm ref}$),
and $\Delta h^{vap}(T_V, P_{\rm ref})$ = heat of vaporization at
($T_V, P_{\rm ref}$). Set
$\Delta H^{vap} := H_o^V - H_o^{LV} \equiv (\rho^V-\rho^L ) h_o^{LV}(T_V) + \rho^V \Delta h^{vap}(T_V, P_{\rm ref})$,
which can be made $> 0$ for any $T_m < T < T_{\rm crit}$
by choosing $h_{\rm ref}$ sufficiently negative, namely
$h_{\rm ref} < -\Delta h^{fus} - c_p^L \, [T_{\rm crit} - T_m]$,
thanks to the fact that $\Delta h^{vap}(T_{\rm crit}) = 0$.
%
\begin{center}
Then $\Lambda^V = (H-H_o^{LV})/\Delta H^{vap}$, $\Lambda^L = 1-\Lambda^V$, $\Lambda^S = 0$,
$\rho = \Lambda^L \rho^L + \Lambda^V \rho^V$,
$H = \rho h = H_o^{LV}(T_V) + \Lambda^V \Delta H^{vap}$,
$h=H/\rho$. \\
\end{center}
\noindent
$\bullet$ \textbf{Vapor phase}\quad ($T > T_V$):\quad
$H > {\bf H_o^V (T_V, P_{\rm ref}) }$.
\begin{center}
Then $\Lambda^V = 1$, $\Lambda^L = 0$, $\Lambda^S = 0$,
$H(T,P) = H_o^{LV}(T_V) + \int_{T_V}^T {C_p^{V}(\tau,P_{\rm ref})} d\tau
+ \int_{P_{\rm ref}}^P { [-T \alpha^V (T, p) +1] }dp$ ,\\
\end{center}
with $\alpha^V$ the thermal expansion coefficient of vapor.
\vspace{2mm}
%\noindent
Since the heat capacity is strictly positive and
$\Delta H^{fus}$ and $\Delta H^{vap}$ are also positive
(by our choice of $h_{\rm ref}$),
the dependence of $H$ on $T$ is monotonic, as shown schematically
in Fig. \ref{fig:3}. Thus $T$ can be found from $H$.
\begin{figure}[ht]
\centering
\includegraphics[width=0.95\textwidth]{fig3} % ,height=3.0in
\label{fig:3} \caption{
Schematic dependence of temperature $T$ and phase fractions
on volumetric enthalpy $H$.}
\end{figure}
The PDE \eqref{eq:3:1} is discretized by a finite volume method,
and explicitly in time.
Briefly, the enthalpy algorithm proceeds as follows:
Knowing the values of $T$ and phase fractions at time step $t_n$,
we evaluate the thermophysical property values, compute the flux
$\vec{Q}$, and update the enthalpy $H$ from the discretized PDE
to time $t_{n+1}$.
From $H$ we find the phase indicators
${\bf{H_o^S}}$, ${\bf{H_o^L}}$, ${\bf{H_o^{LV}}}$, ${\bf{H_o^V}}$,
which determine the temperature and phase fractions as described
above.
\subsection{Knudsen Layer and Vaporization temperature}\label{S:3:2}
%
Whereas melting occurs at a known melt temperature $T_m$
(= 1358 K for copper),
vaporization can occur at any temperature $T_V$, up to the critical
temperature $T_{\rm crit}$ ($\sim 8000$ K for copper),
and must be determined at each time step.
$T_V$ depends on the pressure above the melt, which is coupled to the
Knudsen layer.
Across the Knudsen Layer, treated as a gas dynamic discontinuity
(shock), conservation of mass, momentum and energy impose certain
jump conditions on temperature, pressure, density, and velocity.
There is extensive literature on modeling the Knudsen layer,
at various scales and levels of detail, e.g.
\cite{anisimov68,ytrehus77,knight79,knight81,miotelo95,ohmura98,
ellegaard99,gusarov05,gusarov00}.
A useful treatment is given in Gusarov-Smurov \cite{gusarov05}.
Contrary to most earlier works, their treatment accounts for
evaporation into a background gas. Moreover, it considers condensation
back to the target, occuring when the outer pressure exceeds the
surface pressure.
The KL relations couple the target and plume, and the pressure ratio
across the KL governs what happens at the melt surface.
Let $T_{V}$, $P_{V}$ be the values at the melt surface
(vaporization front),
and $T_{KL}$, $P_{KL}$ denote the values at the outer side of the
Knudsen layer (away from the target, towards the plume).
If the pressure ratio $P_{ratio} := P_{KL}/P_{V} < 1$ then
subsonic evaporation occurs.
If $P_{ratio} > 1$ then condensation occurs
(vaporized material re-condenses onto the melt).
%
Note that in case of ablation into vacuum one can assume sonic
evaporation, in which case the pressure and temperature ratios over
the KL are constant, resulting in the decoupling of the target and
plume processes.
Every quantity experiencing a jump across the Knudsen layer is
ultimately a function of the vaporization temperature $T_{V}$,
which changes in time and needs to be evaluated at each timestep.
We have developed a new approach, based on enthalpy, that includes
the underlying evaporation kinetics.
The KL expressions from \cite{gusarov05}
are applied in an iterative procedure that determines $T_{V}$ in a
vaporizing control volume.
An outline of the approach (too involved to present here in detail):
is as follows:
The updated enthalpy $H$ detects if a control volume is vaporizing
(if ${\bf H_o^V (T_V, P_{\rm ref})} < H < {\bf H_o^V (T_V, P_{\rm ref})}$)
at the current $T_V$.
From the temperature $T_{KL}$, pressure $P_{KL}$, and gas density
$\rho^{KL}$ on the plume side, we find the
recession velocity $v_{V}$ of the vaporization front
(by the Hertz-Knudsen relation and conservation of mass).
From the recession velocity $v_{V}$ we update the vapor fraction of
the (control) volume $Vol$ during a time-step $\Delta t$ via
$\Lambda^V$ = $\Lambda_{old}^V + v_{V} Vol / \Delta t$.
Then we estimate the new enthalpy of the volume as
$H_{new}(T_V) = H_o^{LV}(T_V) + \Lambda^V \Delta H^{vap}(T_V)$,
which should equal $H$. This constitutes an equation
$H_{new}(T_V) = H$ for the unknown $T_V$, which we solve iteratively
(by a bisection type root finder, such as Brent's algorithm).
Beyond the Knudsen layer, the plume of vapor expands rapidly and
ionizes into plasma. The flow is modeled by the Euler equations of gas
dynamics (plus species conservation for electrons, ions, neutrals,
He gas), which must be solved simultaneously with the target
since they are coupled at the Knudsen Layer.
We employ finite volume discretization, explicit time-stepping,
and a central 2nd order high resolution advection scheme.
To save computational time, we use an adaptive grid in the plume
which expands ahead of the plume.
Our numerical scheme allows temperature and pressure dependence
of all thermophysical properties, and vapor properties from
equations of state (such as \cite{TEOS} for copper).
\section{Numerical Simulations}\label{S:4}
We present 1-D simulations of laser ablation of Cu.
The target is initially set at room temperature (300 K).
Above the target we assume a stationary background gas (He), initially
set at atmospheric pressure (1 atm) and room temperature (300 K).
The entire sytem (target+background gas) is irradiated by a laser
at the conditions mentioned in \S\ref{S:3}.
The target, of thickness 12 ${\mu}$m, is discretized along the
axial direction (direction of the laser beam), into a non-uniform
grid of 2000 control volumes, denser near the surface.
In the plume region we use (much) bigger cells,
of length $\Delta z = 5 \times 10^{-7}$ m = 500 nm,
the number of which is adaptively increased to extend the grid
beyond the expanding plume.
The time-step is taken to be $\Delta t = 1 \times 10^{-13}$ s,
respecting the CFL condition in the entire computational domain.
The total simulation time is 30 ns, which corresponds to a single
pulse, cf. Fig. \ref{fig:2}.
Thermodynamic properties of Cu are obtained from the ITTEOS
\cite{TEOS} EOS data (via interpolation).
As a result of the deposited laser energy the target heats up very fast.
The time evolution of the melting and vaporization fronts is shown
in Fig. \ref{fig:4}.
At 7.5 ns melting starts and vaporization is assumed to be minor.
At 10 ns the first cell exceeds the normal boiling point ($T_b=2835K$)
and vaporization becomes sigificant.
Thus the vaporization (surface) temperature $T_V$ and pressure $P_{V}$,
as well as the pressure at the outer side of the KL, $P_{KL}$,
increase quickly as seen in Fig. \ref{fig:6}--\ref{fig:7}.
Since $P_{KL} < P_{V}$ the evaporation is subsonic (see \S\ref{S:3:2}).
\begin{figure}[hpt]
\centering
\includegraphics[width=0.95\textwidth]{fig4}
% , height=4.0in]{fig4}
\label{fig:4}
\caption{ Melt and vaporization depths vs time.}
\end{figure}
When the surface reaches temperatures near the critical state,
the melt becomes transparent. This bleaching of liquid properties
occurs at temperatures about 0.9\,$T_{\rm crit}$
\cite{batanov73}, \cite{karapetyan75}, \cite{yoo2000}.
Here a strong reduction of the electron density in the liquid results
in a drastic change of the optical properties. As the incident laser
energy penetrates through this transparent layer to the underlying
material, the transparency front starts to propagate into the interior
liquid. The transparent surface cell now evaporates at a constant
vaporization temperature $T_V \approx 0.9\,T_{\rm crit}$ as can be seen in
Fig. \ref{fig:5}.
\begin{figure}[hpt]
\centering
\includegraphics[width=0.95\textwidth]{fig5} %,height=3.8in
\label{fig:5} \caption{Vaporization temperature vs time. }
\end{figure}
In this hot dense plume, an amount of electrons is created due to
photoinization of the excited atoms. These electrons absorb additional
energy in, so called, inverse bremssthralung processes
\cite{zelldovich66}.
Due to the aborption of laser energy, additional electrons and ions are
formed resulting in an avalanche effect that keeps on increasing the
plume temperatures to a few 10000 K (here 30000 K): a plasma is created.
At a certain instant, here at $\approx 16$ ns,
the plasma is so dense that it completely shields the target.
As a result, the vaporization temperature and pressure decrease,
and when $P_{V} < P_{KL}$ we arrive at the condensation regime
(Fig. \ref{fig:6}). During this stage, vapor from the plume
recondenses on the target, the recession velocity becomes negative
(Fig. \ref{fig:7}), and the crater depth decreases \eqref{fig:8}.
\begin{figure}[hpt]
\centering
\includegraphics[width=0.95\textwidth]{fig6} % ,height=4.0in
\label{fig:6} \caption{
Vaporization (saturation) pressure $P_{V}$ vs time, and
pressure at the outer side of the Knudsen layer, $P_{KL}$,
for comparison.
}
\end{figure}
\begin{figure}[ht]
\centering
\includegraphics[width=0.75\textwidth]{fig7} %,height=3.2in
\label{fig:7} \caption{
Vaporization (recession) velocity vs time.
}
\end{figure}
After a while, at $\approx 20$ ns, the plume density above the target
decreases due to advection, the laser light reaches again the target
surface raising the surface temperature and pressure,
and we return to the subsonic evaporation regime.
Accordingly, the laser intensity profile at the target takes a bimodal shape (Fig. \ref{fig:2}).
As the applied laser pulse diminishes toward the end,
the amount of laser energy deposited in the plume decreases, the surface
temperature and pressure decrease as well, and material recondenses on
the surface. In the final stage evaporation and condensation compete,
resulting in a recession rate of $\approx 0$ and constant evaporation
depth of $\approx 150$ nm, as seen in Fig. \ref{fig:7} -- \ref{fig:8}.
\begin{figure}[ht]
\centering
\includegraphics[width=0.75\textwidth]{fig8} % ,height=3in
\label{fig:8} \caption{
Vaporization depth vs time. Maximum depth reaches 200 nm shortly after 15 ns (when laser intensity peaks).
}
\vspace{-2mm}
\end{figure}
\section{Conclusion}\label{S:5}
A seamless enthalpy formulation of heat transfer and phase changes
during nanosecond laser ablation of a target is presented.
The formulation is based on Equations of State of the form
$H = H(T, P, phase)$
consistent with the thermochemistry of the material.
It allows full temperature (and/or pressure) dependence of
thermophysical and of optical properties, and can use
available EOS data, up to critical temperature.
Thus, it can realistically describe all the phases of real materials.
The enthalpy model tracks both the solid/liquid and liquid/vapor
fronts and accounts for the kinetic relations at the liquid/vapor
front (Knudsen layer) using an iterative procedure.
It couples naturally to plasma dynamics models in the plume
and the coupled model posseses all the desirable features
listed in \S\ref{S:2}.
Thus, the entire laser ablation process can be simulated,
from heating of the target, to plasma formation and laser absorption,
to the shielding effect of plasma on the target and recondensation.
As far as we know, this is the first model that connects the
thermodynamics and underlying kinetics of this challenging phase
change problem in a self-consistent manner.
Being dimension-independent, it is applicable in 1D, 2D, and 3D
modeling of laser ablation.
As proof of principle, we presented a one-dimensional numerical
implementation for laser ablation of a copper target in helium
background gas at atmospheric pressure, coupled to a gas dynamics
model in the plume (where electons, ions, and He gas are tracked).
As we saw in \S\ref{S:4},
the expected physical behavior of the system is captured very well.
\subsection*{Acknowledgements}
We thank Drs. PR Levashov and KV Khishchenko for providing EOS data
for Cu. We also thank C.C. Garcia of ISAS, Dortmund, for the
drawing of the LA-ICP-MS-setup (Fig. \ref{fig:1}).
Support of the second author by the
Flemish Fund for Scientific Research (FWO Vlaanderen)
is gratefully acknowledged.
\begin{thebibliography}{00}
\bibitem{book} Vasilios Alexiades and Alan D. Solomon,
\emph{Mathematical Modeling of Melting and Freezing Processes},
Hemisphere Publ.Co., 1993.
\bibitem{anisimov68} S. I. Anismov,
"Vaporization of metal absorbing laser radiation",
\emph{Sov. Phys., JETP}, \textbf{27}: 182-183, 1968.
\bibitem{batanov73} V. A. Batanov, F. V. Bunkin, A. M. Prokhorov,
and V. B. Fedorov;
``Evaporation of metals caused by intense optical radiation'',
\emph{Sov. Phys., JETP}, \textbf{36}: 311-322, 1973.
\bibitem{becker02} J. S. Becker,
``Applications of inductively coupled plasma mass spectrometry and laser
ablation inductively coupled plasma mass spectrometry in materials science",
\emph{Spectrochim. Acta, Part B}, \textbf{57}: 1805--1820, 2002.
\bibitem{bleiner05} D. Bleiner,
``Mathematical modelling of laser-induced particulate formation in direct solid microanalysis",
\emph{Spectrochim. Acta Part B}, \textbf{60}(1): 49-64, 2005.
\bibitem{bleiner06} D. Bleiner and A. Bogaerts,
``Computer simulations of laser ablation sample introduction for
plasma-source elemental microanalysis",
\emph{J. Anal. At. Spectrom.}, \textbf{21}: 1161--1174, 2006.
\bibitem{bogaerts03} Annemie Bogaerts, Zhaoyang Chen, Renaat Gijbels, Akos Vertes,
``Laser ablation for analytical sampling: what can we learn from modeling?",
\emph{Spectrochimica Acta B}, \textbf{58}: 1867-1893, 2003.
\bibitem{broekaert02} J. A. C. Broekaert,
\emph{Analytical Atomic Spectrometry with flames and plasmas},
Wiley VCH, 2002.
\bibitem{chen05} Zhaoyang Chen, Annemie Bogaerts;
``Laser ablation of Cu and plume expansion into 1 atm ambient gas",
\emph{J. Appl. Phys}, \textbf{97}: 063305, 2005.
\bibitem{chichkov96} B. N. Chichkov, C. Momma, S. Nolte, F. von Alvensleben,
A. T\"{u}nnermann,
``Femtosecond, picosecond and nanosecond laser ablation of solids",
\emph{Appl. Phys. A}, \textbf{63}: 109--115, 1996.
\bibitem{chrisey94} D. B. Chrisey and G. K. Hubler,
\emph{Pulsed Laser Deposition of Thin Films},
Wiley, 1994.
\bibitem{ellegaard99}
O. Ellegaard, J. Schou, H. M. Urbassek,
``Monte Carlo description of gas flow from laser-evaporated silver",
\emph{Applied Physics A}, \textbf{69}: S577S581, 1999.
\bibitem{garrelie98} F. Garrelie, J. Aubreton, and A. Catherinot,
``Monte Carlo simulation of the laser-induced plasma plume
expansion under vacuum: Comparison with experiments",
\emph{J. Appl. Phys.}, \textbf{83}, 5075, 1998.
\bibitem{gunther02} D. G\"{u}nther,
``Laser ablation-inductively coupled plasma mass spectrometry trends",
\emph{Anal. Bioanal. Chem.}, \textbf{372}: 31--32, 2002.
\bibitem{gusarov00} A. V. Gusarov, A. G. Gnedovets, and I. Smurov,
``Two-dimensional gas-dynamic model of laser ablation in an ambient gas,"
\emph{Applied Surface Science}, \textbf{154-155}: 66-72, 2000.
\bibitem{gusarov05} A. V. Gusarov, I. Smurov,
``Thermal model of nanosecond pulsed laser ablation: Analysis of energy and mass transfer",
\emph{J. Applied Physics}, \textbf{97}: 014307, 2005.
\bibitem{harrington67} R. E. Harrington,
``Application of the Theory of Heat Conduction to the Absorption of
Blackbody Radiation",
\emph{J. Applied Physics}, \textbf{38}: 3266, 1967.
\bibitem{itina97}
T. E. Itina, V. N. Tokarev, W. Marine, and M. Autric,
``Monte Carlo simulation study of the effects of nonequilibrium chemical
reactions during pulsed laser desorption",
\emph{J. Chem. Phys.}, \textbf{106}, 8905, 1997.
\bibitem{itina02} T. E. Itina, J. Hermann, Ph. Delaporte, and M. Sentis,
``Laser-generated plasma plume expansion: Combined continuous-microscopic
modeling", \emph{Phys. Rev. E}, \textbf{66}, 066406, 2002.
\bibitem{itina03} T. E. Itina, J. Hermann, Ph. Delaporte, and M. Sentis,
``Combined continuous--microscopic modeling of laser plume expansion",
\emph{Appl. Surf. Sci.}, \textbf{27}: 208--209, 2003.
\bibitem{karapetyan75} R. V. Karapetyan and A. A. Samokhin,
``Influence of an increase in the transparency on the intense evaporation
of metals by optical radiation",
\emph{Sov. J. Quantum Electron.}, \textbf{4}: 1141--1142, 1975.
\bibitem{knight79} C. J. Knight,
``Theoretical modeling of rapid surface vaporization with back pressure",
\emph{ AIAA Journal}, \textbf{17}(5): 519--523, 1979.
\bibitem{knight81} C. J. Knight,
``Transient Vaporization from a Surface into Vacuum",
\emph{AIAA Journal}, \textbf{20}(7): 950-954, 1981.
\bibitem{kools93} J. C. S. Kools,
``Monte Carlo simulations of the transport of laser-ablated atoms in a diluted gas",
\emph{J. Appl. Phys.}, \textbf{74}, 6401, 1993.
\bibitem{TEOS} P. R. Levashov and K. V. Khishchenko,
ITTEOS 5.8 software for calculation of EOS for metals, 2007.
\bibitem{lindner09} H. Lindner, D. Autrique, C. C. Garcia, K. Niemax
and A. Bogaerts,
``Optimized transport setup for high repetition rate pulse-separated
analysis in laser ablation-inductively coupled plasma mass spectrometry",
\emph{Anal. Chem.}, \textbf{81}: 4241--4248, 2009.
\bibitem{miotelo95} A. Miotello, A. Peterlongo, R. Kelly,
``Laser-pulse sputtering of aluminium: gas-dynamic effects with
recondensation and reflection conditions at the Knudsen layer,"
\emph{Nuclear Instruments and Methods in Physics Research B},
\textbf{101}: 148-155, 1995.
\bibitem{noorbatcha87}
I. NoorBatcha, R. R. Lucchese, and Y. Zeiri,
``Monte Carlo simulations of gas-phase collisions in rapid desorption
of molecules from surfaces",
\emph{J. Chem. Phys.}, \textbf{86}, 5816, 1987.
\bibitem{ohmura98} E. Ohmura, I. Fukumoto, I. Miyamoto,
``Molecular dynamics simulation of laser ablation of metal and silicon",
\emph{Int. J. Jpn. Soc. Prec. Eng.}, \textbf{32}: 248-253, 1998.
Int. J. Jpn. Soc. Prec. Eng. 32
\bibitem{ready71} J. F. Ready,
\emph{Effects of High Power Laser Radiation},
Academic Press, 1971.
\bibitem{vertes93} A. Vertes, R. Gijbels, F. Adams,
\textbf{Laser Ionization Mass Analysis},
Wiley VCH, 1993.
\bibitem{vonAllmen87} M. von Allmen,
\textbf{Laser-Beam Interactions with Materials},
Springer, 1987.
\bibitem{yoo2000} J. H. Yoo, S. H. Jeong, X. L. Mao, R. Greif
and R. E. Russo,
``Evidence for phase-explosion and generation of large particles
during high power nanosecond laser ablation of silicon",
\emph{Appl. Phys. Lett.}, \textbf{76}: 783-785, 2000.
\bibitem{ytrehus77} T. Ytrehus,
``Theory and experiments on gas kinetics in evaporation",
pp.1197-1212 in
\emph{Rarefied Gas Dynamics}, edited by Potter, AIAA, 1977.
\bibitem{zelldovich66}
B. Zelldovich and Yu P. Raizer,
\emph{Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena },
Academic Press, 1966.
\end{thebibliography}
\end{document}