Second-Order Dual Phase Lag Equation. Modeling of Melting and Resolidification of Thin Metal Film Subjected to A Laser Pulse

: The process of partial melting and resolidification of a thin metal film subjected to a high-power laser beam is considered. The mathematical model of the process is based on the second-order dual phase lag equation (DPLE). Until now, this equation has not been used for the modeling of phase changes associated with heating and cooling of thin metal films and the considerations regarding this issue are the most important part of the article. In the basic energy equation, the internal heat sources associated with the laser action and the evolution of phase change latent heat are taken into account. Thermal processes in the domain of pure metal (chromium) are analyzed and it is assumed that the evolution of latent heat occurs at a certain interval of temperature to which the solidification point was conventionally extended. This approach allows one to introduce the continuous function corresponding to the volumetric fraction of solid or liquid state at the neighborhood of the point considered, which significantly simplifies the phase changes modeling. At the stage of numerical computations, the authorial program based on the implicit scheme of the finite difference method (FDM) was used. In the final part of the paper, the examples of numerical computations (including the results of simulations for different laser intensities and different characteristic times of laser pulse) are presented and the conclusions are formulated.


Introduction
Heat transfer through thin films subjected to an ultrafast laser pulse is of vital importance in microtechnology applications and is a reason that the problems related to the fast heating of solids have become a very active research area. The problems of melting/resolidification processes modeling, which may be the result of heating with a laser beam, are also important from the technical point of view. So far, the method using the equation based on the second order model with two delay times for the phase changes modeling was not presented in literature. This was the most important motivation for the authors to undertake research in this area. In Section 5, the comparison of the results obtained with the similar solution on the basis of the first-order dual phase lag equation (DPLE) is presented.
The mathematical model of macroscale heat conduction is based on the parabolic Fourier equation. This equation was formulated under the assumption of instantaneous propagation of the thermal wave in the domain considered. It is obvious that this assumption is not correct, but for the problems concerning the analysis of macroscale heat conduction processes, the obtained results are fully satisfactory. Despite this, the attempts have been made to modify the Fourier equation to a form that better reproduces the real conditions of heat conduction in solids. Thus, about seventy years ago, Cattaneo proposed a modification of the Fourier equation now called the Cattaneo-Vernotte equation. This is the hyperbolic partial differential equation (PDE) and contains the parameter τq called the relaxation time (the lag time of the heat flux in relation to the temperature gradient) [1,2]. Especially important differences between the Fourier model and the real course of thermal processes appear in the case of microscale heat transfer problems. For example, the very high heating rates accompanying the heating of thin metal films with a laser beam mean that the inclusion of the finite value of thermal wave velocity must be taken into account. The deviations appear mainly when the mean free path of the heat carriers becomes comparable to the characteristic length of the domain considered and the time scale of interest becomes comparable to or smaller than the relaxation time of the heat carriers [3,4].
For the analysis of this type of process, the model with two delay times called a dual-phase lag model is presently applied. In addition to the relaxation time, the thermalization time is introduced. The relaxation time τq takes into account the small-scale response in time, while the thermalization time τT takes into account the small-scale response in space [3][4][5][6]. The dual phase lag equation (DPLE) results from the generalized form of the Fourier law where q is a heat flux vector; ∇T is a temperature gradient; λ is a thermal conductivity; and x and t denote the geometrical co-ordinates and time, respectively. Both sides of the last dependence are developed into a power series and, finally (depending on the number of components), the first-or second-order DPLE can be obtained. The literature on equations with two delay times is very extensive (especially in the case of the first order equations), and here we quote only a few important articles. The first publications concerning the model with two delay times appeared in the early nineties of the last century. There may be mentioned, for example, the papers [7][8][9]. Currently, one can already find books devoted to this type of non-Fourier heat conduction model, for example, [10][11][12][13].
In this brief literature review, the selected papers on analytical and numerical solutions of first-order DPLE will be listed. First of all, the works containing the analytical solutions of first-order DPLE (usually 1D tasks) will be mentioned [14][15][16][17][18]. The paper [14] concerns the laser heating of ultra-thin metal film; in the papers [15,16], the bio-heat transfer problems are discussed; while in the paper [18], the multi-layered metal domain is considered.
A much larger number of articles concern the application of numerical methods in the tasks based on the models with two delay times. It should be pointed out that, first of all, different variants of finite difference method (FDM) are applied (see, for example, [19][20][21][22][23][24]). In the paper [19], the numerical model of heating of the double-layered thin film has been applied for the analysis of the thermal deformation process. In the paper [20], the 3D FDM numerical model of the thin metal film heating has been presented. In [21], the explicit scheme of the FDM has been used and the problem of biological tissue freezing process has been discussed. The stability problem of the algorithm of this type is analyzed in [22]. In [23], the problem based on DPLE has been solved using the alternating direction implicit FDM scheme. In the paper [24], the adaptation of typical boundary conditions for non-Fourier equations has been shown. The FDM numerical solutions of the inverse problems can also be found (e.g., [25]).
Literature on the second-order DPLE is not as extensive as for the first-order equations. As an example, the papers [33][34][35][36][37] can be mentioned. The main subject of these papers (except [37]) is related to the construction of algorithms for numerical modeling of problems described by second-order DPLE (the different variants of FDM are used). Additionally, the transformed second-order equations are shown in [20,34,37] and the changed forms are more convenient at the stage of numerical modeling.
At the stage of melting/resolidification modeling, the concept presented in [38] (for the macroscale problems) is applied. The capacity of the internal heat source related to the phase change is proportional to the melting/solidification rate. To define this function (in particular, the volumetric fraction of liquid state fL(T)) in the form of a continuous one, the melting point corresponding to the temperature Tm is conventionally replaced by a certain interval [Tm − ∆T, Tm + ∆T], and then the function discussed can be described by a broken line. For this interval, the substitute thermal capacity is defined and the one domain approach can be used. It should be pointed out that the testing computations show a little impact of the interval ∆T width (within reasonable limits) on the results of numerical simulations.
At the beginning of the part of the article devoted to own research, the assumed form of the dual-phase lag equation and the mathematical formulas determining the laser action and the evolution of phase change latent heat are presented. Both phenomena are taken into account by an introduction to the energy equation of the functions determining the efficiency of internal heat sources. Next, the numerical algorithm based on the implicit scheme of FDM is discussed. In the final part of the paper, the results of numerical computations concerning the heating/cooling process of the thin metal film made of chrome are shown. The conclusions resulting from the performed research are also formulated.

Governing Equations
The starting point for the formulation of the energy equation with delays is the generalized Fourier law (1). To obtain the DPLE, the left and right sides of Equation (1) Let us apply the well-known diffusion equation, namely, where c is a volumetric specific heat and Q(x, t) is a capacity of volumetric internal heat sources. When the components containing the second derivatives (Equation (2)) are taken into account, after mathematical manipulations, Equation (3) takes the form When the melting and resolidification problem is considered, the internal heat source in Equation (4) must contain the term controlling the phase change process. This appropriate source function Qm (x, t) can be defined as (5) where L is the volumetric latent heat phase change and fL(x, t) is the volumetric molten state fraction at the neighborhood of the point considered. The last equation is a generalization of what is well known in the thermal theory of foundry processes definition of Qm (e.g., [38]). The function fL(x, t) is equal to zero at the beginning of the heating process until T1 = Tm − ∆T and fL(x, t) = 1 for T2 > Tm + ∆T (Tm is the melting point). In the interval [Tm − ∆T, Tm + ∆T], the function fL(x, t) changes from 0 to 1 in a linear way (such an assumption is fully acceptable). Generally speaking, the volumetric liquid state fraction is given in the form of broken line, this means The derivative of fL(x, t) with respect to the temperature is equal to 0 for T(x, t) < Tm − ∆T and T(x, t) > Tm + ∆T, while between the border temperatures dfL(x, t)/dT = 1/2∆T. Thus, the source term Qm(x, t) acts only for T(x, t) from the interval [Tm − ∆T, Tm + ∆T], and then Let us introduce the piece-vise constant function C(T) (8) where c1 and c2 are the volumetric specific heats of the solid and liquid states, respectively. Then, Equation (4) can be written as follows Thermal conductivity λ in Equation (9) is defined just like the parameter C(T). The mathematical formula determining the intensity of the internal heat source Ql (x, t) resulting from the laser action can be taken in the form [39] where I0 [J/m 2 ] is the laser intensity, tp [s] is the characteristic time of laser pulse, δ [m] is the optical penetration depth, R is the reflectivity of the irradiated surface, rD [m] is the laser beam radius, and x3 is a vertical axis. The derivatives of Ql with respect to time can be found analytically.
On the outer surface of the system, the adiabatic conditions are assumed (the external heat flux is taken into account in the appropriate source function). The mathematical form of the Neumann boundary condition for the second-order DPLE is as follows [36] where n is a normal outward vector and (in the case considered) qb (x, t) = 0, of course.
The mathematical model is also supplemented by the initial conditions where Tp is an initial temperature.

Mathematical Description of 1D Problem
At the stage of numerical modeling, the 1D problem was considered and the basis for the construction of the FDM algorithm is the following system of equations and conditions: − energy equation for thin metal film domain − source function Ql adiabatic boundary conditions

Let
( , ) For the transition t f−1 → t f (f ≥ 3), the approximate form of Equation (13) resulting from the introduction of the assumed differential quotients is as follows while s = f, s = f − 1, or s = f − 2. Thus, 2  1  1  2  1  2  3   2  3   2  2  1  1  1  1  1  1  1  1  2  1 2   2  1  1  1  1  1  1  1  1  2  1 2   τ  2  3  3  τ  2   2  2τ  τ  λ  λ  λ  λ  4   2τ 2τ and after mathematical manipulations, one has 1 1 Different temporal-spatial meshes were considered; see Table 1. For each combination of mesh steps, the temperature after 80 ps on the heated surface was recorded. As can be seen, the result does not depend much on the time step, while the number of nodes is important. Analysis of the results obtained showed that the satisfactory accuracy provides the geometrical mesh containing n = 1000 nodes (h = 0.1 nm) and time step ∆t = 0.0005 ps (see Table 1, column 5).  Figure 1 shows the temperature courses at the points x = 0 (heated surface) and x = 20 nm. The results were compared to those published in [24], where the first order DPL equation and slightly different partial melting model were used. As one can see, the differences are small, but visible. For the second-order model (solid lines), the maximum temperature of the heated surface occurs after 28.8 ps and is equal to 2323.28 K, while for the first-order model (dashed lines), the maximum temperature occurs after 28 ps and equals 2307.20 K. In Figure 2, the fragments of heating/cooling curves for different widths of ΔT are shown. One can see that the results are similar. For ΔT = 3 K, the maximum temperature is the highest; for ΔT = 7 K, the maximum temperature is the lowest; while the biggest differences between them are of the order 10 K. After the resolidification process, the cooling curves for all cases are very similar. Further computations were carried out for ΔT = 3 K. In Figure 3, the temperature courses on the irradiated surface for the laser intensity I0 = 3825 J/m 2 and different characteristic times of laser pulse tp are presented. The impact of changes in this parameter on the course of the heating/cooling curves is clearly visible. It is obvious that an increase in intensity of the laser pulse causes a more rapid heating process. Interesting, however, are the effects of changes of characteristic time tp. Shortening this time increases the heating rate of the metal film and maxima of the cooling/heating curves shift to the left. At the same time, the maximum values of temperatures increase slightly. A detailed analysis of Equation (10) determining the time-dependent efficiency of the laser-related internal heat source shows that such an effect could be expected. Finally, Figure 4 shows the changes of molten layer thicknesses in time. The successive curves correspond to tp = 2, 6, and 10 ps. For tp = 2 ps, the maximal depth is equal to g = 17.6 nm; for tp = 6 ps, it is equal to g = 15.3 nm; while for tp = 10 ps, it is equal to g = 13 nm. Changes in the characteristic time of laser pulse cause a similar effect to that seen in the previous figure. The depth of the molten layer increases with the reduction of time tp and the process is more dynamic. Such analysis and information obtained on its basis can often be useful in engineering practice.

Conclusions
The mathematical model, the numerical algorithm, and exemplary results of the computations concerning the heating of thin metal film subjected to the laser beam are discussed, wherein the problem is described using the second-order DPLE. The laser beam power is so high that, in the domain under consideration, the partial melting and then resolidification of the material take place. To our knowledge, the application of the second-order DPLE to solve this type of the problem has not yet been presented. The task is treated as a non-linear one and the thermophysical parameters of the material are temperature-dependent. The melting point is substituted by a certain interval of temperature. This assumption allows one to introduce the continuous function determining the local and temporary volumetric liquid state fraction of the metal. The influence of the width of this interval on the results of numerical calculations is also examined. The laser action and the evolution of latent heat are taken into account by the introduction of two internal heat sources to the energy equation.
At the stage of numerical modeling, the 1D problem was considered, but the generalization of the algorithm and computer program implementing the numerical computations for the larger number of dimensions is not difficult. The computer program is based on the implicit scheme of the finite difference method. The algorithm is unconditionally stable and the transition from time t to t + ∆t requires the solution of the system of linear equations, whose main matrix is a three-diagonal one. In this place, the very effective Thomas algorithm was used.
The obtained results are in line with the expectations; their comparison with the solution of a similar problem described by the first-order equation shows visible, but slight differences.
We are planning the following further research in this scope: numerical algorithm and computer program for 3D problems; -numerical algorithm and computer program for axially-symmetrical tasks (this geometry is very convenient because of the typical shape of the function describing the laser action); -adaptation of the algorithm presented in this work for modeling the ablation process; -research on other approaches to phase changes modeling.