First and second order dual phase lag equation. Numerical solution using the explicit and implicit schemes of the finite difference method

In the paper the different variants of the dual phase lag equation (DPLE) are considered. As one knows, the mathematical form of DPLE results from the generalization of the Fourier law in which two delay times are introduced, namely the relaxation time τq and the thermalization one τT. Depending on the order of development of the left and right hand sides of the generalized Fourier law into the Taylor series one can obtain the different forms of the DPLE. It is also possible to consider the others forms of equation discussed resulting from the introduction of the new variable or variables (substitution). In the paper a thin metal film subjected to a laser pulse is considered (the 1D problem). Theoretical considerations are illustrated by the examples of numerical computations. The discussion of the results obtained is also presented.


Introduction
The macroscopic heat conduction model (the Fourier model) results from the assumption of instantaneous propagation of the thermal wave in the domain under considerations. It is known that this approach is not admittedly correct, but in the most of problems related to the heat conduction in the macro scale is fully satisfactory.
Seventy years ago Cattaneo [1] formulated an equation in which the delay time (relaxation time τq) of the heat flux in relation to the temperature gradient was taken into account. This hyperbolic PDE is known as the Cattaneo-Vernotte equation. In the case of the macroscopic problems the Cattaneo-Vernotte equation is used as a model of heat transfer processes for the nontypical materials with a complex internal structure (e.g. biological tissue) [2]).
The very high heating rates typical for the microscale heat transfer cause that the inclusion of the finite value of thermal wave velocity must be somehow taken into account.
The very popular model describing the microscale heat transfer is based on the dual phase lag equation (DPLE) -e.g. [3,4]. The final form of DPLE results from the generalized form of the Fourier law in which both the relaxation time and the thermalization one are introduced (e.g. [5]). Depending on the number of terms in the Taylor series expansion of this law, the different forms of the dual phase lag equation can be obtained, in particular, the first and the second order DPLE and also the mixed variants. Some simple problems described by the DPLE and supplemented by the appropriate boundary and initial conditions can be solved analytically. As an example the paper [6] can be mentioned. The solution concerns a thin metal film heating caused by a laser action, while this action is taken into account by the introduction of the artificial internal heat source.
The numerous approximate solutions of the different problems described by the first order DPLE can be found in literature, e.g. [7][8][9]. The numerical solutions concerning the second order DPLE (based on the finite difference method) are the subject of works [10][11][12][13][14]. In the papers [12,13] the mixed variants of DPLE are also considered.
Below, two variants of DPLE modified forms will be discussed. They are the result of a substitutions causing the appearance in the energy equation of the fragment which is identical as in the classical Fourier equation. The others components can be treated as the internal volumetric heat source. The computer program based on the explicit scheme of finite difference method is not complicated, additionally the stability condition has a simple form. At the stage of numerical modeling the 1D problem is considered. It results from the thermal aspects of the task discussed.

Governing equations
The starting point for the formulation of the dual phase lag equation (DPLE) is the generalized Fourier law in the form Dual phase lag equation must be supplemented by the boundary conditions. Used in this paper the Neumann condition takes a the form (6) or (the second order model) where qb(x, t) is a given boundary heat flux, the sign "+" corresponds to x = 0 (left boundary), while sign "-" corresponds to x = G (right boundary), G is the metal film thickness.
The initial conditions are in the form Additionally (for equations (4) and (5) where Tp is an initial temperature, u(x) and v(x) are the known functions.

Modifications of DPLE
In the paper [8] a certain modification of the first order DPLE using the substitution technique is presented. In this paper we consider the second order DPLE (as in [12]) and the energy equation is extended by the introduction of the internal heat sources. Let from which results that Introducing (11) and (12) into (4) gives where a = λ/c is a thermal diffusion coefficient. This form of DPLE equation is interesting because the first two components (containing the function U) correspond to the classical diffusion equation, while the remaining ones can be treated as an internal heat source (thanks to which the formulation of the stability criterion for the explicit differential scheme used here is very simple). The boundary conditions should be also modified. In particular The last formula corresponds to the Neumann condition, of course. The initial condition takes a form (c.f. equations (8) and (9)) 2 τ 0 : Thus, at first the equation (13) supplemented by boundary conditions (13) and initial condition (15) should be solved. Next, using the formula (10), the temporary temperature field can be determined.
Another approach involves the introduction of two substitutions, namely [13] Introducing formula (16) into (17) one obtains the dependence (10) and then equation (4) can be written in the fom The last formula can be differentiated with respect to time and then or (c.f. (17) By introducing (19) into (21) one obtains As in the case of equation (13), the first two components containing the function U correspond to the classical diffusion equation, while the remaining ones can be treated as an artificial internal heat source.
The equation (23) is supplemented by the modified boundary conditions and the initial condition (15). After solving the problem (23), (24), (15), the ordinary differential equation (17) with the initial condition (c.f. formula (16)) 2 τ 0 : is solved. Finally, the solution of the ordinary differential equation (16) with the initial condition T(x,0) = Tp, allows one to determine the temporary temperature field.

Finite Difference Method
The problems presented previously have been solved using the different variants of the finite difference method (explicit and implicit schemes [12,16]).
All above presented variants of DPLE have been modeled using the in-house computer programs based on the FDM, while two selected solutions will be shown and compared here. In particular, the thermal processes proceeding in the domain of thin metal film subjected to the laser pulse are considered. The first solution applies equation (23) with the internal volumetric heat source resulting from the laser action. The energy equation is supplemented by the Neumann conditions (in the adiabatic version) and the appropriate initial conditions. At the stage of numerical modeling the explicit scheme of the FDM has been used. One can see, that equation (23) has a form which is convenient for numerical modeling. It contains the simple differential operators without the mixed derivatives. The same physical problem has been solved using the equation (4) and the adequate boundary-initial conditions. In this version the implicit scheme of the FDM has been applied.
The following approximation of equation (23) is proposed: The well known stability criterion for parabolic equations determines the critical time step.
The FDM form of the boundary conditions is the following (x = 0): and (x = G): and where I is the laser intensity, tp is the characteristic time of laser pulse, δ is the optical penetration depth, R is the reflectivity of the irradiated surface and β = 4 ln2. For x = 0 and x = G the no-flux boundary conditions are assumed: qb(0, t) = qb(G, t)=0.
In Figure 1 the comparison of results obtained using the algorithms discussed is presented. Solid lines correspond to the explicit scheme of FDM for the equation (23) with the boundary conditions (24) and initial condition (15), while the symbols correspond to the implicit scheme of FDM for the equation (4) with the boundary conditions (7) and initial conditions (8), (9). In this Figure the temperature histories at the points x = 0 (irradiated surface), x = 20 nm and x = 40 nm are shown. As can be seen, the results are almost the same. For example, for the first and the second variants of FDM, the maximum temperatures at the irradiated surface are equal to 306.536 K and 306.535 K, respectively.   Figure 3 illustrates the heating/cooling curves at the irradiated surface for a shorter laser pulse tp = 0.05 ps and laser intensities 13.7 J/m 2 and 50 J/m 2 , respectively. For I = 13.7 J/m 2 the maximum temperatures are equal to 315.12 K and 308.28 K, respectively, while for I = 50 J/m 2 the maximum temperatures are equal to 355.17 K and 330.22 K. As previously, for the higher laser intensity the differences between the solutions obtained using the first and second order DPL are greater. For the other intensities: I = 100 J/m 2 : maximum temperatures: 410.34 K and 360.44 K; I = 200 J/m 2 : maximum temperatures: 520.68 K and 420.88 K; I = 300 J/m 2 : maximum temperatures: 631.02 K and 481.33 K.
It should be noted that the results presented in Figures 2 and 3 have been obtained using the implicit scheme of the FDM.

Conclusions
In the paper the second order dual phase lag models are considered. The results of numerical computations show that obtained solutions are practically the same, despite that the final forms of the DPLE are different. The models resulting from the introduction of certain substitutions are more complicated than the classical one, but on the other hand, the algorithms using the FDM are essentially simpler. The results obtained have been compared with the solutions of the first order DPLE and the differences are clearly visible (see: previous chapter).