**4. Numerical Model Based on FDM**

Let *T<sup>f</sup> <sup>i</sup>* = *T*(*xi*, *f*Δ*t*), where *xi* = *ih*, *<sup>i</sup>* <sup>=</sup> 0, 1, ... , *<sup>n</sup>* and *<sup>f</sup>* <sup>=</sup> 0, 1, ... , *<sup>F</sup>*. Taking into account the initial conditions (17), one has

$$\begin{array}{lcl} T\_i^0 = T\_P\\ \frac{T\_i^1 - T\_i^0}{\Delta t} = \frac{1}{\varepsilon\_1} Q\_l(\mathbf{x}\_i; \mathbf{0}) \rightarrow \ T\_i^1 = T\_i^0 + \frac{\Delta t}{c\_1} Q\_l(\mathbf{x}\_i; \mathbf{0})\\ \frac{T\_i^2 - 2T\_i^1 + T\_i^0}{\left(\Delta t\right)^2} = \frac{1}{\varepsilon\_1} \left. \frac{\partial Q\_l(\mathbf{x}; t)}{\partial t} \right|\_{t=0} \rightarrow \ T\_i^2 = 2T\_i^1 - T\_i^0 + \frac{\left(\Delta t\right)^2}{c\_1} \left. \frac{\partial Q\_j(\mathbf{x}\_i; t)}{\partial t} \right|\_{t=0} \end{array} \tag{18}$$

For the transition *t <sup>f</sup>*−<sup>1</sup> <sup>→</sup> *t <sup>f</sup>* (*f* <sup>≥</sup> 3), the approximate form of Equation (13) resulting from the introduction of the assumed differential quotients is as follows

$$\begin{split} &C\_{i}^{f-1} \left[ \frac{T\_{i}^{f} - T\_{i}^{f-1}}{\Delta t} + \pi\_{q} \frac{T\_{i}^{f} - 2T\_{i}^{f-1} + T\_{i}^{f-2}}{\left(\Delta t\right)^{2}} + \frac{\pi\_{q}^{2}}{2} \frac{T\_{i}^{f} - 3T\_{i}^{f-1} + 3T\_{i}^{f-2} - T\_{i}^{f-3}}{\left(\Delta t\right)^{3}} \right] = \\ &\left[ \frac{\partial}{\partial x} \left(\lambda \frac{\partial T}{\partial x}\right) \right]\_{i}^{f} + \frac{\pi\_{T}}{\Delta t} \left\{ \left[ \frac{\partial}{\partial x} \left(\lambda \frac{\partial T}{\partial x}\right) \right]\_{i}^{f} - \left[ \frac{\partial}{\partial x} \left(\lambda \frac{\partial T}{\partial x}\right) \right]\_{i}^{f-1} \right\} + \\ &\frac{w\_{T} \pi\_{T}^{2}}{2 \left(\Delta t\right)^{2}} \left\{ \left[ \frac{\partial}{\partial x} \left(\lambda \frac{\partial}{\partial x} \frac{T}{\partial x}\right) \right]\_{i}^{f} - 2 \left[ \frac{\partial}{\partial x} \left(\lambda \frac{\partial}{\partial x} \frac{T}{\partial x}\right) \right]\_{i}^{f-1} + \left[ \frac{\partial}{\partial x} \left(\lambda \frac{\partial}{\partial x} \frac{T}{\partial x}\right) \right]\_{i}^{f-2} \right\} + Z\_{i}^{f} \end{split} \tag{19}$$

where

$$\begin{aligned} \left[\frac{\partial}{\partial x}\left(\lambda \frac{\partial^{\top} T}{\partial \dot{x}}\right)\right]\_{\dot{\mathbf{i}}}^{\mathcal{S}} &= \frac{1}{\hbar} \Big[ \left(\lambda \frac{\partial}{\partial \dot{x}} \frac{T}{\lambda}\right)\_{\dot{\mathbf{i}} + 0.5}^{\mathcal{S}} - \left(\lambda \frac{\partial^{\top} T}{\partial \dot{x}}\right)\_{\dot{\mathbf{i}} - 0.5}^{\mathcal{S}} \right] = \frac{1}{\hbar} \Big(\lambda\_{\dot{i} + 0.5}^{\mathcal{S}} \frac{T\_{i + 1}^{\ast} - T\_{i}^{\ast}}{\hbar} - \lambda\_{\dot{i} - 0.5}^{\mathcal{S}} \frac{T\_{i}^{\ast} - T\_{i - 1}^{\ast}}{\hbar} \Big) = \\ \frac{1}{\hbar} \Big(\frac{\lambda\_{i}^{\ast} + \lambda\_{i + 1}^{\ast}}{2} \frac{T\_{i + 1}^{\ast} - T\_{i}^{\ast}}{\hbar} - \frac{\lambda\_{i + 1}^{\ast} + \lambda\_{i}^{\ast}}{2} \frac{T\_{i}^{\ast} - T\_{i - 1}^{\ast}}{\hbar} \Big) = \frac{1}{2\hbar^{2}} \Big[ \left(\lambda\_{i}^{\ast} + \lambda\_{i + 1}^{\ast}\right) \Big(T\_{i + 1}^{\ast} - T\_{i}^{\ast}\Big) + \left(\lambda\_{i - 1}^{\ast} + \lambda\_{i}^{\ast}\right) \Big(T\_{i - 1}^{\ast} - T\_{i}^{\ast}\Big) \Big] \end{aligned} \tag{20}$$

while *s* = *f*, *s* = *f* − 1, or *s* = *f* − 2.

Thus,

*Tf <sup>i</sup>* <sup>−</sup>*Tf*−<sup>1</sup> *i* <sup>Δ</sup> *<sup>t</sup>* + τ*<sup>q</sup> Tf <sup>i</sup>* <sup>−</sup>2*Tf*−<sup>1</sup> *<sup>i</sup>* <sup>+</sup>*Tf*−<sup>2</sup> *i* (Δ*t*) <sup>2</sup> <sup>+</sup> <sup>τ</sup><sup>2</sup> *q* 2 *Tf <sup>i</sup>* <sup>−</sup>3*Tf*−<sup>1</sup> *<sup>i</sup>* <sup>+</sup>3*Tf*−<sup>2</sup> *<sup>i</sup>* <sup>−</sup>*Tf*−<sup>3</sup> *i* (Δ*t*) <sup>3</sup> = 2(Δ*t*) 2 +2τ*T*Δ *t*+τ<sup>2</sup> *T* 4*Cf*−<sup>1</sup> *<sup>i</sup> <sup>h</sup>*2(Δ*t*) 2 λ*f*−<sup>1</sup> *<sup>i</sup>* <sup>+</sup> <sup>λ</sup>*f*−<sup>1</sup> *i*+1 -*T f <sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>T</sup> <sup>f</sup> i* + - λ*f*−<sup>1</sup> *<sup>i</sup>*−<sup>1</sup> <sup>+</sup> <sup>λ</sup>*f*−<sup>1</sup> *i* -*T f <sup>i</sup>*−<sup>1</sup> <sup>−</sup> *<sup>T</sup> <sup>f</sup> i* − 2τ*T*Δ *t*+2τ<sup>2</sup> *T* 4*Cf*−<sup>1</sup> *<sup>i</sup> <sup>h</sup>*2(Δ*t*) 2 λ*f*−<sup>1</sup> *<sup>i</sup>* <sup>+</sup> <sup>λ</sup>*f*−<sup>1</sup> *i*+1 -*T <sup>f</sup>*−<sup>1</sup> *<sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>T</sup> <sup>f</sup>*−<sup>1</sup> *i* + - λ*f*−<sup>1</sup> *<sup>i</sup>*−<sup>1</sup> <sup>+</sup> <sup>λ</sup>*f*−<sup>1</sup> *i* -*T <sup>f</sup>*−<sup>1</sup> *<sup>i</sup>*−<sup>1</sup> <sup>−</sup> *<sup>T</sup> <sup>f</sup>*−<sup>1</sup> *i* + τ2 *T* 4*Cf*−<sup>1</sup> *<sup>i</sup> <sup>h</sup>*2(Δ*t*) 2 λ*f*−<sup>2</sup> *<sup>i</sup>* <sup>+</sup> <sup>λ</sup>*f*−<sup>2</sup> *i*+1 -*T <sup>f</sup>*−<sup>2</sup> *<sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>T</sup> <sup>f</sup>*−<sup>2</sup> *i* + - λ*f*−<sup>2</sup> *<sup>i</sup>*−<sup>1</sup> <sup>+</sup> <sup>λ</sup>*f*−<sup>2</sup> *i* -*T <sup>f</sup>*−<sup>2</sup> *<sup>i</sup>*−<sup>1</sup> <sup>−</sup> *<sup>T</sup> <sup>f</sup>*−<sup>2</sup> *i* <sup>+</sup> *<sup>Z</sup><sup>f</sup> i Cf*−<sup>1</sup> *i* (21)

and after mathematical manipulations, one has

$$A\_i^f T\_{i-1}^f + B\_i^f T\_i^f + C\_i^f T\_{i+1}^f = D\_i^f \tag{22}$$

where

*Af <sup>i</sup>* = − [2(Δ*t*) 2 +2τ*T*Δ *t*+τ<sup>2</sup> *T*] λ*f*−<sup>1</sup> *<sup>i</sup>*−<sup>1</sup> <sup>+</sup>λ*f*−<sup>1</sup> *i* 4*Cf*−<sup>1</sup> *<sup>i</sup> <sup>h</sup>*2(Δ*t*) <sup>2</sup> , *Bf <sup>i</sup>* <sup>=</sup> <sup>2</sup>(Δ*t*) 2 +2τ*q*Δ *t*+τ<sup>2</sup> *q* 2(Δ*t*) <sup>3</sup> + [2(Δ*t*) 2 +2τ*T*Δ *t*+τ<sup>2</sup> *T*] λ*f*−<sup>1</sup> *<sup>i</sup>*−<sup>1</sup> <sup>+</sup>2λ*f*−<sup>1</sup> *<sup>i</sup>* <sup>+</sup>λ*f*−<sup>1</sup> *i*+1 4*Cf*−<sup>1</sup> *<sup>i</sup> <sup>h</sup>*2(Δ*t*) <sup>2</sup> , *Cf <sup>i</sup>* = − [2(Δ*t*) 2 +2τ*T*Δ *t*+τ<sup>2</sup> *T*] λ*f*−<sup>1</sup> *<sup>i</sup>*+<sup>1</sup> <sup>+</sup>λ*f*−<sup>1</sup> *i* 4*Cf*−<sup>1</sup> *<sup>i</sup> <sup>h</sup>*2(Δ*t*) <sup>2</sup> , *Df <sup>i</sup>* <sup>=</sup> <sup>2</sup>(Δ*t*) 2 +4τ*<sup>q</sup>* Δ*t*+3τ<sup>2</sup> *q* 2(Δ*t*) <sup>3</sup> *<sup>T</sup>f*−<sup>1</sup> *<sup>i</sup>* <sup>−</sup> <sup>2</sup>τ*<sup>q</sup>* <sup>Δ</sup>*t*+3τ<sup>2</sup> *q* 2(Δ*t*) <sup>3</sup> *<sup>T</sup>f*−<sup>2</sup> *<sup>i</sup>* <sup>+</sup> <sup>τ</sup><sup>2</sup> *q* 2(Δ*t*) <sup>3</sup> *<sup>T</sup>f*−<sup>3</sup> *<sup>i</sup>* − 2τ*T*Δ*t*+2τ<sup>2</sup> *T* 4*Cf*−<sup>1</sup> *<sup>i</sup> <sup>h</sup>*2(Δ*t*) 2 λ*f*−<sup>1</sup> *<sup>i</sup>* <sup>+</sup> <sup>λ</sup>*f*−<sup>1</sup> *i*+1 -*T <sup>f</sup>*−<sup>1</sup> *<sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>T</sup> <sup>f</sup>*−<sup>1</sup> *i* + - λ*f*−<sup>1</sup> *<sup>i</sup>*−<sup>1</sup> <sup>+</sup> <sup>λ</sup>*f*−<sup>1</sup> *i* -*T <sup>f</sup>*−<sup>1</sup> *<sup>i</sup>*−<sup>1</sup> <sup>−</sup> *<sup>T</sup> <sup>f</sup>*−<sup>1</sup> *i* + τ2 *T* 4*Cf*−<sup>1</sup> *<sup>i</sup> <sup>h</sup>*2(Δ*t*) 2 λ*f*−<sup>2</sup> *<sup>i</sup>* <sup>+</sup> <sup>λ</sup>*f*−<sup>2</sup> *i*+1 -*T <sup>f</sup>*−<sup>2</sup> *<sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>T</sup> <sup>f</sup>*−<sup>2</sup> *i* + - λ*f*−<sup>2</sup> *<sup>i</sup>*−<sup>1</sup> <sup>+</sup> <sup>λ</sup>*f*−<sup>2</sup> *i* -*T <sup>f</sup>*−<sup>2</sup> *<sup>i</sup>*−<sup>1</sup> <sup>−</sup> *<sup>T</sup> <sup>f</sup>*−<sup>2</sup> *i* <sup>+</sup> *<sup>Z</sup><sup>f</sup> i Cf*−<sup>1</sup> *i* (23)

The approximation of boundary conditions (16) is as follows

$$\cdot \quad \text{ for } \ge 0;$$

$$\frac{T\_1^f - T\_0^f}{h} + \frac{\pi\_T}{\Delta t} \left( \frac{T\_1^f - T\_0^f}{h} - \frac{T\_1^{f-1} - T\_0^{f-1}}{h} \right) + \frac{\pi\_T^2}{2(\Delta t)^2} \left( \frac{T\_1^f - T\_0^f}{h} - 2\frac{T\_1^{f-1} - T\_0^{f-1}}{h} + \frac{T\_1^{f-2} - T\_0^{f-2}}{h} \right) = 0 \tag{24}$$


$$\frac{T\_n^f - T\_{n-1}^f}{h} + \frac{\tau\_T}{\Delta t} \left( \frac{T\_n^f - T\_{n-1}^f}{h} - \frac{T\_n^{f-1} - T\_{n-1}^{f-1}}{h} \right) + \frac{\tau\_T^2}{2(\Lambda t)^2} \left( \frac{T\_n^f - T\_{n-1}^f}{h} - 2\frac{T\_n^{f-1} - T\_{n-1}^{f-1}}{h} + \frac{T\_n^{f-2} - T\_{n-1}^{f-2}}{h} \right) = 0 \tag{25}$$

or

$$\begin{cases} -\left[2\left(\Delta t\right)^2 + 2\tau\_T \Delta t + \tau\_T^2\right] T\_0^f + \left[2\left(\Delta t\right)^2 + 2\tau\_T \Delta t + \tau\_T^2\right] T\_1^f = \\ \left(2\tau\_T \Delta t + 2\tau\_T^2\right) \left(T\_1^{f-1} - T\_0^{f-1}\right) - \tau\_T^2 \left(T\_1^{f-2} - T\_0^{f-2}\right) \end{cases} \tag{26}$$

and

$$\begin{cases} -\left[2\left(\Delta t\right)^2 + 2\tau\_T \Delta t + \tau\_T^2\right] T\_{n-1}^f + \left[2\left(\Delta t\right)^2 + 2\tau\_T \Delta t + \tau\_T^2\right] T\_n^f = \\ \left(2\tau\_T \Delta t + 2\tau\_T^2\right) \left(T\_n^{f-1} - T\_{n-1}^{f-1}\right) - \tau\_T^2 \left(T\_n^{f-2} - T\_{n-1}^{f-2}\right) \end{cases} \tag{27}$$

The transition from *t <sup>f</sup>* <sup>−</sup> <sup>1</sup> to *t <sup>f</sup>* (*f* <sup>≥</sup> 3) requires the solution of the system of Equations (22), (26), and (27) with a three-band main matrix that can be solved quickest using the Thomas algorithm [40].

#### **5. Results of Computations**

The computations were performed for the thin metal film of the thickness *G* = 100 nm made of chromium. The surface *x* = 0 is subjected to the laser pulse with the parameters *R* = 0.93, *I*<sup>0</sup> = 3825 J/m2, δ = 15.3 nm, and *tp* = 10 ps—c.f. Formula (15). The following values of chromium parameters are assumed: thermal conductivities λ<sup>1</sup> = 93 W/(m K), λ*<sup>2</sup>* = 35 W/(m K), volumetric specific heats *c*<sup>1</sup> = 3.2148 MJ/(m3 K), *c*<sup>2</sup> = 2.79276 MJ/(m<sup>3</sup> K) [24], relaxation time τ*<sup>q</sup>* = 0.136 ps, thermalization time τ*<sup>T</sup>* = 7.86 ps, melting temperature *Tm* = 2180 K, and volumetric heat of fusion *L* = 2904 MJ/m3.

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

**Table 1.** Temperature of the irradiated surface after 80 ps for different time steps and number of nodes.


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.

**Figure 1.** Temperature histories at the points *x* = 0 and *x* = 20 nm; comparison of the results obtained using the first-order (dashed line [24]) and second-order dual phase lag (DPL) (solid line) for *I*<sup>0</sup> = 3825 J/m2, *tp* = 10 ps.

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.

**Figure 2.** Comparison of obtained temperature courses for different intervals Δ*T* (on the surface *x* = 0), (1) Δ*T* = 3 K, (2) Δ*T* = 5 K, (3) Δ*T* = 7 K.

In Figure 3, the temperature courses on the irradiated surface for the laser intensity *I*<sup>0</sup> = 3825 J/m<sup>2</sup> 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.

**Figure 3.** Heating/cooling curves on the surface *x* = 0 for different characteristic times of laser pulse. *tp*: (1) *tp* = 2 ps, (2) *tp* = 6 ps, (3) *tp* = 10 ps (*I*<sup>0</sup> = 4500 J/m2).

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.

**Figure 4.** Changes in the thickness of the molten layer for different characteristic times of laser pulse, (1) *tp* = 2 ps, (2) *tp* = 6 ps, (3) *tp* = 10 ps.

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.

## **6. 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:



**Author Contributions:** Conceptualization, E.M. and B.M.; formal analysis, E.M. and B.M.; software, E.M.; investigation, E.M. and B.M., writing—original draft, B.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.
