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.

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 Mazumder3 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).

FIG. 1.

Surface temperature instability during pulsed-laser ablation and its significance for ablation of deeper craters. The color bars in the image plots denote the cell volume or fill fraction in a fixed coordinate that ranges from 0 to 1. Rapid variations in the surface temperature are mainly caused by the cell volume.

FIG. 1.

Surface temperature instability during pulsed-laser ablation and its significance for ablation of deeper craters. The color bars in the image plots denote the cell volume or fill fraction in a fixed coordinate that ranges from 0 to 1. Rapid variations in the surface temperature are mainly caused by the cell volume.

Close modal

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.

FIG. 2.

(a) Normalized (extrapolated) surface temperature stability with cell index number across an ablated crater surface and (b) extrapolated surface temperature correlation with fill fraction, particularly at very low values of fill fraction.

FIG. 2.

(a) Normalized (extrapolated) surface temperature stability with cell index number across an ablated crater surface and (b) extrapolated surface temperature correlation with fill fraction, particularly at very low values of fill fraction.

Close modal

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.

The heat conduction equation implemented originally2 can be written as
H t = z ( κ T z ) + r ( κ T r ) + Q ( t , z , r ) ,
(1a)
H = T 0 T ρ ( u ) C p ( u ) d u ,
(1b)
where 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 Cp(T) are the density and specific heat capacity, respectively. Both ρ and Cp 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.6 
The boundary condition of an energy or heat flux is prescribed on the surface, i.e., at the edges of the surface grid cells z = Z s ( r , t ), which moves as ablation occurs,
( κ T z , κ T r ) z = Z s ( r , t ) = ( F z , F r ) ,
(2)
where F z and F r are the surface heat fluxes in z and r directions, respectively. The heat flux at the crater surface includes cooling rates induced by the radiation and evaporation processes, respectively,
κ T z | z = Z s ( r , t ) = σ s b ( ε s T s 4 ε a m b T a m b 4 ) w L l v ,
(3)
where σ s b is the Stefan–Boltzmann constant, L l v is the latent heat associated with the liquid-vapor phase change, and ( ε s , ε a m b ) and ( T s , T a m b ) are the emissivity and temperature for the surface and ambient air, respectively.
For the purposes of transparently illustrating the occurrence of the said surface temperature instability, we employ a simple (controlled) 1D treatment and add analytical formulations necessary for simulating ablation here. The surface recession speed (w) is related to the metal evaporation rate, and we adopt the Hertz–Knudsen equation for the evaporation rate ( φ ) of a melting metal. This leads to a linear relation between w, φ, and the surface pressure (ps),
w = m φ ρ l = α e m 2 π k B T s p s ρ l ,
(4)
where k B is the Boltzmann constant, m is the atomic mass, α e (=0.95) is the coefficient of evaporation, and ρ l is the liquid metal density. Furthermore, p s and surface temperature ( T s ) generally are related by the following empirical relation:
log ( p s / atm ) = A + B T s 1 + C log T s ,
(5)
where the constants (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 Ts(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).

FIG. 3.

(a) Surface temperature and (b) surface cell index/location, vs time.

FIG. 3.

(a) Surface temperature and (b) surface cell index/location, vs time.

Close modal
FIG. 4.

(a) Surface pressure and (b) surface recession (ablation) velocity, vs time.

FIG. 4.

(a) Surface pressure and (b) surface recession (ablation) velocity, vs time.

Close modal
In this illustrative example, we assume an aluminum substrate and prescribe the input energy flux F z ( t ) used in Eq. (8) for a pulsed-laser beam at 10.6 μm in the form
F z ( t ) = R s F 0 exp [ ( t 1.5 τ ) 2 / τ 2 ] ,
(6)
with F 0 = 2 × 10 12 W m 2, τ = 5 μs, and the surface reflectivity R S = 0.926. The corresponding coefficient of evaporation, α e, used in Eq. (4) is set to 0.95. We set a varying spatial step size in the z-direction Δ z j + 1 = 1.03 Δ z j with Δ z 1 = 0.5 μ m and j = 1 , 2 , , 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.
An effective way to resolve the “underfill” problem is to introduce an ablation-front-tracking frame of reference by employing a moving-coordinate transform [Eqs. (7) and (8) and Fig. 5],
t = t ,
(7a)
r = r ,
(7b)
z = z w t = z w ( t , r ) t ,
(7c)
where the surface receding velocity w ( t , r ) = w ( t , r ) is related to the other surface physical parameters by Eq. (4). Because the new coordinate ( t , z , r ) follows the surface receding velocity, the cell size of the surface cell in ( t , z , r ) will remain unchanged, which effectively solves the underfill problem in a fixed coordinate for which the size of the surface cell will gradually decrease to zero until the moving surface hits an underlying new cell. Correspondingly, the heat conduction Eq. (1a) can be transformed from the fixed coordinate ( t , z , r ) to a moving coordinate ( t , z , r ),
H t = w H z + ( t w t ) H z + z ( κ T z ) + r ( κ T r ) + Q ( t , z , r ) r ( κ t w r T z ) t w r z ( κ ( T r t w r T z ) ) .
(8)
Most importantly, the boundary conditions (2) and (3) are greatly simplified and can now be specified on a fixed boundary of a horizontal surface ( z = 0 ),
κ T z | z = 0 = F z ( t ) σ s b ( ε s T s 4 ε a m b T a m b 4 ) w L l v .
(9)
FIG. 5.

Schematic diagram of a fixed vs a moving coordinate in the z-direction at different time steps. The size of the surface cell in a fixed coordinate changes with time due to an ablation process. On the other hand, the cell size remains unchanged in a moving coordinate. Note that the positive z-direction is downward.

FIG. 5.

Schematic diagram of a fixed vs a moving coordinate in the z-direction at different time steps. The size of the surface cell in a fixed coordinate changes with time due to an ablation process. On the other hand, the cell size remains unchanged in a moving coordinate. Note that the positive z-direction is downward.

Close modal

Because the boundary condition (9) has been prescribed on a fixed surface of z = 0, we can adopt an operator splitting method7 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.

Comparing Eq. (8) with Eq. (1a), we note that the moving-coordinate transform (7) introduces four new terms in the original 2D heat conduction equation. In this paper, we retain the first two terms on the right-hand side of the transformed Eq. (8) because our numerical experiments indicate that these two terms are essential. However, the last two of the four terms may be neglected for the present 1D analysis because these two terms are of a higher order for representing rz correlations. It may be noted that other heat conduction algorithms with a receding surface4,8 retain only the first term [ w ( H / z ) ] 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
H t = w H z + ( w ( log t ) ) H z + z ( κ T z ) + Q ( t , z , r ) .
(10)
The heat conduction, Eq. (10) [or Eqs. (1) and (8)], can be solved by the Dufort–Frankel pseudo-implicit scheme that is stable for constant coefficients.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 transition6 
H = { T 0 T ρ ( u ) C p ( u ) d u , T < T m T 0 T m ρ C p d u + T m T ρ C p d u + L s l , T > T m ,
(11)
where Lsl is the enthalpy jump across the melt point Tm associated with the solid–liquid phase change. The above discontinuity in H–T relation also leads to an undetermined value H at Tm. We note that, in general, other model parameters such as the heat conduction coefficient κ and (ρ, Cp) are also discontinuous functions of temperature across the melt front defined by Tm, which also leads to the nonuniqueness of their values at Tm. Physically, the above discontinuous H–T relation (11) at Tm further leads to a corresponding discontinuity in temperature variation with time.
Mathematically, the discontinuities in physical parameters across an internal surface mean that the governing heat conduction equation needs to be solved in two domains separately with a set of compatibility conditions that ensures matching solutions on both sides of the internal surface. One example of using such an approach to solve the heat conduction equation including laser-induced ablation was given by Chan and Mazumder3 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 Tm. Specifically, we rewrite the enthalpy function (11) in the following analytic form for any T:
H = T 0 T ρ ( u ) C p ( u ) d u + H ( T ) ,
(12)
with H*(T) expressed in the following form:
H ( T ) = L s l 2 [ 1 + tanh ( T T m Δ ) ] L s l U H ( T T m ) , as Δ 0 ,
(13)
where U H ( x ) denotes a Heaviside unit step function and L s l = L m ρ m is the enthalpy jump at the melt point due to the latent heat ( L m ) and the discontinuous changes in metal density ( ρ m ). The so-called Moving Frame Solver is a combination of the moving-coordinate transform (7) that leads to a new heat conduction Eq. (8) or Eq. (10) and the introduction of a sigmoid function (13) to deal with the parameter jumps in the system. The specification Δ depends on how well one likes to resolve the structure of a melt front. In practice, given a grid size of a numerical model and for a temperature difference among neighboring cells, say, Δ T 10 K, one may set Δ = Δ T / 100 = 0.1 K to resolve the melt front fine structure. Physically, Δ plays the role of the solid fraction in the previous treatment of the melting front.1,6 One advantage of using Eqs. (12) and (13) to replace Eq. (11) is that the term dH/dT that is required in solving Eq. (8) or Eq. (10) by using an implicit scheme can also be expressed as an analytic function
d H d T = [ ρ ( T ) C ( T ) ] + d H d T Δ H Δ T + L s l 2 Δ sec h 2 ( T T m Δ ) .
(14)
The first term on the right-side of Eq. (14) represents the slow variation of 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
1 2 Δ sec h 2 ( T T m Δ ) δ ( T T m ) , as Δ 0 ,
(15)
where δ ( x ) is the Dirac–δ 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 Tm.
Note that H / t = ( d H / d T ) ( T / t ); as a result, the heat conduction Eq. (10) can be rewritten as the one that only includes one dependent variable T,
T t ( w + w log t ) T z ( d H d T ) 1 z ( κ T z ) = ( d H d T ) 1 Q ( t , z , r ) .
(16)

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 dH/dT in the third term is given by Eq. (14) that approximately equals ( ρ C p ) 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.

FIG. 6.

Validation of the implicit thermal solver and sigmoid function by comparison with analytical solutions: (a) laser-induced thermal profile without melting (Ref. 9) and (b) laser-induced thermal profile with melting and melt front motion (Stefan problem) (Ref. 11).

FIG. 6.

Validation of the implicit thermal solver and sigmoid function by comparison with analytical solutions: (a) laser-induced thermal profile without melting (Ref. 9) and (b) laser-induced thermal profile with melting and melt front motion (Stefan problem) (Ref. 11).

Close modal

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 = 0. A direct comparison of ablation surface parameters demonstrates the effectiveness of the Moving Frame model in mitigating the instability issue.

FIG. 7.

Comparisons of the numerically simulated surface parameters as a function of time between a fixed-coordinate frame approach (blue curves) and a moving-coordinate frame approach (red curves): (a) surface temperature, (b) surface pressure, (c) surface recession velocity, and (d) ablation crater depth.

FIG. 7.

Comparisons of the numerically simulated surface parameters as a function of time between a fixed-coordinate frame approach (blue curves) and a moving-coordinate frame approach (red curves): (a) surface temperature, (b) surface pressure, (c) surface recession velocity, and (d) ablation crater depth.

Close modal

The laser absorption depth is d ∼ 10 nm,6 and the spatial grid size, Δ z 1 = 0.5 μ 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 Ts (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 Tm = 933 K for aluminum) moving toward the interior of the substrate. Specifically, a melt front characterized by a near-constant temperature at Tm moves from a surface cell (cell-1) to the underlying cell-2 around time t = 1.6 μs in a jumping fashion in a discretized numerical scheme.

FIG. 8.

Surface temperature Ts and the temperatures at five cells, i.e., cell center temperatures, close to the surface in the first 3 μm seconds for a base run using the moving-frame coordinate scheme. A melt front with a near-constant temperature at Tm moves from surface cell (cell-1) to the neighboring internal cell (cell-2).

FIG. 8.

Surface temperature Ts and the temperatures at five cells, i.e., cell center temperatures, close to the surface in the first 3 μm seconds for a base run using the moving-frame coordinate scheme. A melt front with a near-constant temperature at Tm moves from surface cell (cell-1) to the neighboring internal cell (cell-2).

Close modal

Figure 9 shows the crater depth (δc) and the melt front depth (δsl) 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.

FIG. 9.

Crater depth (δ) (blue line) and melt front depth (δsl) (red line) evolution.

FIG. 9.

Crater depth (δ) (blue line) and melt front depth (δsl) (red line) evolution.

Close modal

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 × 10 12 W m 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.

FIG. 10.

(a) Surface temperature and (b) melt front depth as a function of radial distance and time after an integration of 5 μs forced by a temporal-spatial Gaussian pulse.

FIG. 10.

(a) Surface temperature and (b) melt front depth as a function of radial distance and time after an integration of 5 μs forced by a temporal-spatial Gaussian pulse.

Close modal

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.

The authors thank Anthony Mark for assistance with generating Figs. 1 and 2.

The authors have no conflicts to disclose.

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).

1.
C. P.
Grigoropoulos
,
T. D.
Bennett
,
J.-R.
Ho
,
X.
Xu
, and
X.
Zhang
, “
Heat and mass transfer in pulsed-laser-induced phase transformations
,”
Adv. Heat Transfer
28
,
75
144
(
1996
).
2.
C.
Sawyer
,
K.
Iyer
,
X.
Zhu
,
M.
Kelly
,
D.
Luke
, and
D.
Amdahl
, “
Two-dimensional laser-induced thermal ablation modeling with integrated melt flow and vapor dynamics
,”
J. Laser Appl.
29
,
022212
(
2017
).
3.
C. L.
Chan
and
J.
Mazumder
, “
One-dimensional steady-state model for damage by vaporization and liquid expulsion due to laser-material interaction
,”
J. Appl. Phys.
62
,
4579
4586
(
1987
).
4.
Y.
Wang
,
N.
Shen
,
G. K.
Befekad
,
L.
Crystal
, and
C. L.
Pasiliao
, “
Modeling pulsed laser ablation of aluminum with finite element analysis considering material moving front
,”
Int. J. Heat Mass Transfer
113
,
1246
1253
(
2017
).
5.
B.
Wang
,
Y.
Huang
,
J.
Jiao
,
H.
Wang
,
J.
Wang
,
W.
Zhang
, and
L.
Sheng
, “
Numerical simulation on pulsed laser ablation of the single-crystal superalloy considering material moving front and effect of comprehensive heat dissipation
,”
Micromachines
12
,
225
(
2021
).
6.
C. P.
Grigoropoulos
,
Transport in Laser Microfabrication: Fundamentals and Applications
(
Cambridge University. Press
,
2009
),
400
p.
7.
W. H.
Press
,
S. A.
Teukolsky
,
W. T.
Vetterling
, and
B. P.
Flannery
,
Numerical Recipes in Fortran. The Arts of Scientific Computing
2nd ed. (
Cambridge University Press
,
1992
),
963
p.
8.
N. M.
Bulgakova
and
A. V.
Bulgakov
, “
Pulsed laser ablation of solids: Transition from normal vaporization to phase explosion
,”
Appl. Phys. A
73
,
199
208
(
2001
).
9.
G. J.
Haltiner
and
R. T.
Williams
,
Numerical Prediction and Dynamic Meteorology
, 2nd ed. (
John Wiley & Sons
,
1980
),
477
p.
10.
J. H.
Bechtel
, “
Heating of solid targets with laser pulses
,”
J. Appl. Phys.
46
,
1585
1593
(
1975
).
11.
H. S.
Carslaw
and
J. C.
Jaeger
,
Conduction of Heat in Solids
, 2nd ed. (
Oxford University Press
,
1959
),
510
p.