An accepted approach to computing laser-induced peak surface temperature is to employ the enthalpy formulation of the transient heat conduction equation [Grigoropoulos *et al.*, Adv. Heat Transfer **28**, 75–144 (1996); Sawyer *et al.*, J. Laser Appl. **29**, 022212 (2017)]. This approach is generally implemented using an explicit numerical scheme to solve the thermal transport equation. While it offers the advantage of modeling the solid-melt phase transition automatically, the approach results in instability-like behavior in the computed surface temperature. When laser-induced ablation becomes significant, the heating rate in the surface cell becomes unrealistically large. This results in spikes in the computed peak surface temperature due to large errors in calculating the heating rate. In this paper, we present a new approach, which we refer to as the Moving Frame Solver, that employs a moving-coordinate frame of reference, located at the receding evaporating surface. We also use an analytical representation for the phase transition region of the enthalpy-temperature relationship. The Moving Frame Solver combined with an implicit scheme leads to a stable solution without surface temperature, pressure, or velocity spikes. In other words, any instability in these computed parameters due to use of an explicit scheme (such as Dufort–Frankel) has been eliminated. Details of the new thermal solver and example calculations are presented. Numerical experiments suggest that the surface cell size needs to be small, ∼0.1 *μ*m, to obtain a highly accurate solution with a typical metal such as aluminum. Using the Moving Frame Solver with a refined grid near the surface, but coarse elsewhere, enables accurate and stable surface temperature computation.

## I. INTRODUCTION

Predictive material response modeling of pulsed-laser ablation is necessary for understanding the sensitivity of various parameters controlling ablation efficiency. In recent years, we have been developing a two-dimensional (2D), axisymmetric volume of fluid (VOF)-based computational framework, DEIVI (directed energy illumination and visualization) for predicting pulsed-laser ablation of metals.^{1,2} The DEIVI code broadly follows the approaches of Grigoropoulos *et al*.^{1} and Chan and Mazumder^{3} and directly models several metal material response phenomena in a fully coupled manner: the melt transition, motion of the solid-melt interface, Navier–Stokes surface melt flow, melt evaporation, evaporative gasdynamics of the vapor using both a linear Lax–Wendroff and piecewise parabolic Colella–Woodward scheme, and plasma transition of the vapor and laser absorption by plasma. While the previous approaches were one-dimensional (1D) in their treatment of the metal substrate and simplified the treatment of the melt response or ignored it entirely, we have implemented all of these material response aspects in a complete 2D framework. Also, neither of the previous approaches modeled plasma generation and absorption.

While the DEIVI code has shown excellent predictive capability for ablation involving relatively shallow crater depths in aluminum (and for short duration nanosecond pulses), the code sometimes shows a tendency for exhibiting surface temperature instability when simulating illumination profiles that generate deeper craters [>O(10 *μ*m)] (Fig. 1).

In the area of pulsed-laser ablation modeling, it is a known and longstanding challenge that actual surface temperatures during ablation (surface recession) are not known and cannot be measured easily. Extreme transient values, e.g., 4000–6000 K in aluminum, can be attained over nanoseconds and involve rapid material phase change. Consequently, there is no available means of validating surface temperature computations directly. Models are forced to rely on computing the substrate (solid phase) thermal field and use some type of extrapolation procedure to compute the highly transient, nonequilibrium temperature at the ablating surface as it is critical for driving all aspects of ablation via evaporation and recoil-pressure-driven melt expulsion. Despite the criticality of surface temperature computation, there is little discussion in the literature about addressing this longstanding challenge for pulsed-laser ablation modeling and complications encountered in implementation of the many models published in the open literature. We seek to address this gap and present a mitigation approach with this paper.

In the case of the DEIVI code, the root cause of the issue noted above was ultimately traced to a type of instability that arises in the surface temperature computation when the receding surface approaches the bottom edge of a cell in a fixed grid (Fig. 2). This situation corresponds to a state where the subject VOF cell has a very low (near zero) fill fraction.

Initially, we attempted some revisions of the extrapolation scheme used to compute surface temperature. These alternate schemes were unable to mitigate the surface temperature instability issue, which is as expected, since the root cause lies with computing temperature in a highly underfilled VOF cell. One way to deal with this core problem is to ensure that every surface cell always has a high fill fraction. By employing a coordinate frame of reference that is always located at the recessed surface location,^{3} one can ensure that the surface cell is always full (fill fraction = 1). In this paper, we describe our approach to updating the DEIVI thermal solver by employing such a model, which we refer to as the Moving Frame model. In implementing this, we also approximate the solid-melt phase change discontinuity in the enthalpy-temperature relationship by use of a sigmoid function to enable a fully implicit solution of the recast heat equation.

Compared with other recent efforts that employ a coordinate frame of reference that moves with the receding surface during pulsed-laser ablation,^{4,5} the DEIVI code is unique (and more advanced) owing to its direct modeling of a melt transition and a Navier–Stokes melt flow model in a VOF framework. This naturally allows crater shape evolution resulting from ablation not just due to evaporation but also melt expulsion; in a fully coupled manner. In the other models referenced above, the substrate is modeled as a solid until the material reaches a specified ablation temperature (which the authors do not specify) in a finite element (FE) domain. A separate simplified algorithm computes surface shape evolution resulting from an analytical Hertz–Knudsen equation (for evaporation only). While the objective of these other approaches with adopting a Moving Frame approach is to improve accuracy by attempting to eliminate unnecessary laser-induced heating of evaporated material, our objective is to mitigate instability in the surface temperature computation.

Section II A describes the current thermal model and illustrates the instabilities that arise in surface temperature and ablation parameters, consequentially using a simple analytical model. Section II B describes the Moving Frame approach used to mitigate the surface temperature instability. Section II C describes the analytical function employed for the enthalpy-temperature curve, including the discontinuity across the solid-melt phase transition, to enable solution with just a single variable (temperature). Section III presents results showing the mitigation of instabilities using the Moving Frame approach.

## II. MODEL DESCRIPTION

### A. DEIVI thermal model

^{2}can be written as

*t*is the time and

*z*and

*r*denote the depth and radius in a cylindrical coordinate, respectively. We use the convention for the vertical coordinate (

*z*) that it points into the metal substrate, originating at the initial surface. In Eq. (1),

*T*is the temperature, with $ T 0$ being a reference temperature (298 K),

*κ*is the heat conductivity coefficient, and

*ρ*(

*T*) and

*C*(

_{p}*T*) are the density and specific heat capacity, respectively. Both

*ρ*and

*C*are functions of temperature. The term $Q(t,z,r)$ in Eq. (1) represents the heating rate generated by the volumetric absorption of the laser energy in the current problem. Equation (1) has been implemented in the enthalpy form to conveniently track the movement of the solid-melt interface.

_{p}^{6}

*z*and

*r*directions, respectively. The heat flux at the crater surface includes cooling rates induced by the radiation and evaporation processes, respectively,

*w*) is related to the metal evaporation rate, and we adopt the Hertz–Knudsen equation for the evaporation rate $(\phi )$ of a melting metal. This leads to a linear relation between

*w*, $\phi $, and the surface pressure (

*p*),

_{s}*m*is the atomic mass, $ \alpha e$ (=0.95) is the coefficient of evaporation, and $ \rho l$ is the liquid metal density. Furthermore, $ p s$ and surface temperature $( T s)$ generally are related by the following empirical relation:

*A*,

*B*,

*C*) are determined experimentally.

Equations (1)–(5) form a closed system for which the metal temperature distribution *T*(*t*,*z*,*r*) including its surface temperature *T _{s}*(

*t*,

*r*) can be solved for given external heat source

*Q*(

*t*,

*z*,

*r*) and metal property parameters such as

*ρ*(

*T*),

*κ*(T),

*C*(

*T*), and $ L l v$. This approach, which is similar to the one implemented in DEIVI, exhibits surface temperature instability (Figs. 3 and 4). When ablation becomes significant (in relation to surface cell size), the heating rate in the surface cell becomes unrealistically large when the surface cell index undergoes a change (in the fixed coordinate frame of the reference used).

*μ*m in the form

*τ*= 5

*μ*s, and the surface reflectivity $ R S=0.926$. The corresponding coefficient of evaporation, $ \alpha e$, used in Eq. (4) is set to 0.95. We set a varying spatial step size in the

*z*-direction $\Delta z j + 1=1.03\Delta z j$ with $\Delta z 1=0.5\mu m$ and $j=1,2,\u2026,120$ to deal with the fact that the temperature variation in the region near the surface is greater than that away from the surface and requires a higher spatial resolution to resolve its structure.

### B. Moving frame approach

Because the boundary condition (9) has been prescribed on a fixed surface of $ z \u2032=0$, we can adopt an operator splitting method^{7} to solve Eq. (8) for which we only need to develop a numerical scheme for a 1D heat conduction equation and then alternatively apply it to Eq. (8) in *z* and *r* directions. Under such a circumstance, we may set the heat flux at the lateral boundary to zero $ F r=0$.

*r*–

*z*correlations. It may be noted that other heat conduction algorithms with a receding surface

^{4,8}retain only the first term $[ w ( \u2202 H / \u2202 z \u2032 )]$ on the right-hand side in Eq. (8). In summary, the 1D heat conduction equation in the

*z*-direction that we solve in the Moving Frame model is given by

### B. Enthalpy-temperature parameterization for a fully implicit solution

^{9}Our numerical experiments indicate that a more effective and stable numerical scheme to solve the equation with variable coefficients, i.e., Eq. (10), is by the use of a fully implicit numerical scheme because the triangular matrix can be efficiently inverted by a forward and backward substitution method.

^{7}In order to adopt the forward and backward substitution method, however, there should exist only one dependent variable in the equation to be solved directly. However, there are two dependent variables,

*T*and

*H*, in Eq. (9) that are related by Eq. (1b). Further, the functional relation between

*T*and

*H*is not a continuous function due to the latent heat release associated with the melt transition

^{6}

*L*is the enthalpy jump across the melt point

_{sl}*T*associated with the solid–liquid phase change. The above discontinuity in

_{m}*H–T*relation also leads to an undetermined value

*H*at

*T*. We note that, in general, other model parameters such as the heat conduction coefficient

_{m}*κ*and (

*ρ*,

*C*) are also discontinuous functions of temperature across the melt front defined by

_{p}*T*, which also leads to the nonuniqueness of their values at

_{m}*T*. Physically, the above discontinuous

_{m}*H–T*relation (11) at

*T*further leads to a corresponding discontinuity in temperature variation with time.

_{m}^{3}where different moving-coordinates and analytical solutions were obtained and matched in the solid and liquid domains separately. An alternate method to deal with the internal discontinuity surface that is especially effective in numerically solving Eq. (10) [or Eq. (1)] is to approximate an internal discontinuity surface by a rapidly but continuously varying internal surface so that the two domains may be merged into one. In Grigoropoulos

*et al.*,

^{1}the mathematical discontinuity in the

*H–T*relation (11) was eliminated by introducing an additional parameter in the numerical scheme that calculates a solid/liquid fraction within a finite cell volume at the melt front to give a definitive value of

*H*at the cell containing the melt point. In this paper, we introduce a new parameterization scheme that uses a sigmoid function to uniformly represent all discontinuity parameters across the melt temperature

*T*. Specifically, we rewrite the enthalpy function (11) in the following analytic form for any

_{m}*T*:

*H**(

*T*) expressed in the following form:

^{1,6}One advantage of using Eqs. (12) and (13) to replace Eq. (11) is that the term d

*H*/d

*T*that is required in solving Eq. (8) or Eq. (10) by using an implicit scheme can also be expressed as an analytic function

*H*with respect to

*T*and can be evaluated numerically by finite difference. The second term is the parametrization of the discontinuous variation of

*H*with respect to

*T*at the melt point $ T m$, and it has been parameterized as an analytic function. We note that

*δ*function. Thus, our parameterization of using Eq. (14) to replace Eq. (11) is not only simple but also mathematically well-founded. The same sigmoid function can also be applied to the heat conduction coefficient

*κ*, which has a large discontinuity change across

*T*.

_{m}*T*,

Inspecting the left-hand side of the above heat conduction equation, we note that the added second term comes from the introduction of the moving coordinate, whereas the coefficient d*H*/d*T* in the third term is given by Eq. (14) that approximately equals $ ( \rho C p ) \u2212 1$ away from the melt point. At or near the melt front, the second term in Eq. (14) is far greater than the first term. This will greatly reduce the magnitude of the third term in Eq. (16), which leads to a very slow change in temperature and the corresponding melt front position with time as expected.

In summary, our Moving Frame Solver consists of two changes to the widely used fixed-coordinate heat conduction equation solver: (i) a moving-coordinate transform (7) to yield a new Eq. (10) in the moving *z*′-coordinate and (ii) a sigmoid function transform (13) that is applied to all discontinuously varying parameters to merge different physical domains into one. Both these changes were validated by comparison with standard known solutions for time-dependent laser-induced heating and cooling,^{10} as well as melt front motion (Stefan problem),^{11} shown in Fig. 6.

## III. RESULTS

To demonstrate the merits of the two changes individually, we first show numerical results that adopt the first improvement only. In Fig. 7, we show the numerically simulated surface parameters as a function of time (red curves) for the system given by Eqs. (4), (5), and (10) with a boundary condition of Eq. (9) prescribed on the moving surface of $ z \u2032=0$. A direct comparison of ablation surface parameters demonstrates the effectiveness of the Moving Frame model in mitigating the instability issue.

The laser absorption depth is *d* ∼ 10 nm,^{6} and the spatial grid size, $\Delta z 1=0.5\mu m$, is much greater. As a result, the laser heating only occurs at the surface cell and its magnitude per grid of material at the surface cell will be inversely proportional to the grid size as the ablation occurs. If one still adopts an inter-cell heat conduction scheme with a constant cell size, the mismatch between the laser heating and the conductive heat transport toward the internal neighboring cell at the surface cell results in an unrealistically large surface temperature as the surface cell shrinks to zero as shown in Fig. 7(a) by the blue curve. Such a mismatch not only leads to a surface temperature vacillation associated the vacillation of the surface cell size but also leads to a significant overestimate of the average surface temperature. The corresponding surface pressure [Fig. 7(b)] and surface recession velocity [Fig. 7(c)] calculated from Eqs. (4) and (5) also show the spike features caused by the unrealistic surface temperature. When the neighboring cell becomes a new surface cell, the surface temperature solved by using a fixed-coordinate scheme drops but the numerical errors continue to accumulate. While the crater depth evolution [Fig. 7(d)] does not show instability, its value remains inaccurate due to the cumulative errors in the solved recession surface velocity as shown in Fig. 7(c). In the moving-frame coordinate (Fig. 5), the cell size is independent of the surface recession velocity associated with the ablation process. As a result, the surface heating rate by the absorption of laser energy and its conductive transport to the internal neighboring cell are always appropriately balanced, yielding a smooth solution.

When the second improvement is also included, the model becomes capable of tracking the melt front motion also. Figure 8 shows the surface temperature *T _{s}* (which is receding in relation to the initial surface location as ablation occurs) and the temperatures at five cells, i.e., cell center temperatures, immediately beneath the moving surface in the first 3

*μ*m seconds. The figure clearly shows the melt front (with a near-constant temperature at

*T*= 933 K for aluminum) moving toward the interior of the substrate. Specifically, a melt front characterized by a near-constant temperature at

_{m}*T*moves from a surface cell (cell-1) to the underlying cell-2 around time

_{m}*t*= 1.6

*μ*s in a jumping fashion in a discretized numerical scheme.

Figure 9 shows the crater depth (*δ _{c}*) and the melt front depth (

*δ*) as a function of time. The region sandwiched between the blue and red lines is the melt potion of the domain. In reality, the melt front depth will be less due to lateral melt flow and recoil-pressure-driven expulsion, which is not included in the present 1D model.

_{sl}By adopting an operator splitting method to solve the 2D axisymmetric heat conduction equation, the 1D Moving Frame solver can be easily incorporated into a 2D module for solving the heat conduction, Eq. (8). In Fig. 10, we show an example of the 2D plots for the surface temperature and the melt front depth as a function of radial distance and time after an integration of 5 *μ*s. The fields are forced by a temporal-spatial Gaussian pulse defined by Eq. (6) with a peak irradiance of $ F 0=0.4\xd7 10 12 W m \u2212 2$ and a temporal half-width of *τ* = 0.5 *μ*s plus a spatial half-width of 200 *μ*m. The temporal-spatial distributions of the surface temperature and the melt front depth are expected for the given temporal-spatial Gaussian pulse. Again, in practice, the lateral melt flow forced by the recoil pressure drives expulsion that could more effectively affect the melt front depth.

## IV. CONCLUSIONS

We have developed a new algorithm that we refer to as the Moving Frame Solver as a way to mitigate issues with surface temperature calculation related to surface recession in VOF-based computation of laser ablation. This new algorithm introduces a moving-coordinate transform such that laser heating at the receding surface can be calculated more accurately by eliminating extrapolation from interior (cell center) temperature values in low fill fraction (underfilled) VOF cells. Our heat conduction equation in the moving-coordinate frame Eq. (10) contains two extra new terms that are essential and comparable in magnitude. In addition, a sigmoid function representation is employed for the enthalpy transition across the solid-melt temperature transition so that the heat conduction equation in the two domains with different states can be solved uniformly in one composite domain as shown by Eq. (16). Validation of our implicit numerical scheme to solve the heat conduction equation and the sigmoid function representation is demonstrated by comparison with well-known exact solutions. A laser ablation computation example is presented to illustrate the efficacy of the new approach and algorithm.

## ACKNOWLEDGMENT

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Xun Zhu:** Conceptualization (equal); Formal analysis (lead); Methodology (equal); Software (lead); Visualization (lead); Writing – original draft (lead). **Kaushik A. Iyer:** Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). **Darren Luke:** Funding acquisition (equal); Methodology (equal); Resources (equal); Supervision (equal); Visualization (equal); Writing – review & editing (equal).

## REFERENCES

*Transport in Laser Microfabrication: Fundamentals and Applications*

*Numerical Recipes in Fortran. The Arts of Scientific Computing*

*Numerical Prediction and Dynamic Meteorology*

*Conduction of Heat in Solids*