Sensitivity of a Process for Heating Thin Metal Film Described by the Dual-Phase Lag Equation with Temperature-Dependent Thermophysical Parameters to Perturbations of Lag Times

Abstract

In the paper, an equation with two delay times (dual-phase lag Equation (DPLE)) in a version that takes into account the dependence of thermophysical parameters (volumetric specific heat and thermal conductivity) on temperature is considered. In particular, an analysis of the sensitivity of transient temperature field in relation to disturbances in delay times (the relaxation and thermalization times) is performed. The sensitivity model concerns the process of heating an ultrathin metal layer with a laser beam. First, the equation with two delay times in the case of temperature-dependent thermophysical parameters is presented. Next, the sensitivity equations with respect to delay times are derived using the direct method. The algorithms for solving the basic and sensitivity tasks are also briefly presented. At the stage of computations, an authorial program based on the implicit scheme of a finite-difference method is developed. In the final part of the paper, examples of numerical solutions (for layers made from gold and nickel) are presented. The research conducted here shows that disturbances in the temperature field are clearly visible and depend, on the one hand, on the thermophysical parameters of the material, and on the other hand, on the intensity of heating with an external heat source.

1. Introduction

The significant development of various types of microtechnologies in recent years has resulted in the need for an increasingly detailed analysis of thermal processes occurring on the microscale. Research in this field concerned first the formulation of appropriate mathematical models and then the methods for solving the problems using analytical, semi-analytical and, most significantly, numerical approaches.

The article presents the results of research concerning DPLE (the dual-phase lag equation) in the version discussed in [1]. In the cited work, a numerical algorithm for solving a version of the dual-phase lag equation in which thermophysical parameters are temperature-dependent is discussed. Here, we analyze the sensitivity of such a model to delay times.

The variants of DPLE are often used to model microscale heat-transfer problems. It may also be mentioned here that two other frequently used mathematical models for the description of heat transport on the microscale are the two-temperature models (parabolic and hyperbolic [2,3,4]) because under some assumptions, the dual-phase lag equation can be obtained from the parabolic two-temperature model [5]. Additionally, on the basis of this model, the values of delay times for various metal materials can be determined analytically.

The mathematical forms of equations with two delay times result from the adoption of a generalized Fourier law in the form [6,7,8]

$$q\left(X,t+{\tau }_q\right)=-\lambda (T)\mathrm{grad}T\left(X,t+{\tau }_T\right)$$

(1)

where q is the heat flux, T is the temperature, X is the geometrical coordinates, t is time, τ_q (relaxation time) is the delay in the generation of the heat flux and the associated local temperature gradient, τ_T (thermalization time) is the delay in the generation of the temperature gradient caused by the heat conduction through small-scale structures, and λ is the thermal conductivity. In sum, the temperature gradient at point X for time t + τ_T corresponds to the heat flux at the same point for time t + τ_q. The lag times for typical metals are of the order of a few to tens of picoseconds; therefore, in the case of macro-scale problems, these values are irrelevant, while for the heat transfer on a microscale (especially for subdomains subjected to intense ultra-fast external thermal influences), the delay times must be taken into account at the stage of problem formulation. This factor is why this paper undertakes the analysis of the impact of perturbations of these parameters on the space-time temperature fields in the domain considered.

The first equation in which the delay time was introduced is the Cattaneo–Vernotte equation [9]. In this equation, the authors took into account the final velocity of thermal-wave propagation.

The basic version of DPLE results from the expansion of the left and right sides of Equation (1) into a Taylor series with accuracy to the first derivative, e.g., [6,7,8]. In the literature, equations with delay times resulting from the expansion of the generalized Fourier’s law with accuracy to the second-order derivatives [10,11,12], as well as the mixed variants [13], are also considered. In turn, the integro-differential form of DPLE is presented in [14], while the models with lag times in which derivatives with respect to time are of non-integer order are presented, among others, in [15,16]. Mathematical models based on the DPLE are commonly used for the analysis of microscale heat-transfer processes. In the literature, one can find the few works in which the equation with three lag times is considered, e.g., [17,18]. In the papers on delay-time equations, the limits that ensure the physical correctness of the model are also analyzed. In most of them, the mutual relationships between delay times are formulated mathematically [10,19,20].

In recent years, the several papers [1,21,22] concerning the more general form of DPLE for variable, temperature-dependent thermophysical parameters (volumetric specific heat and thermal conductivity) have been presented. In [21], the variability of the thermal conductivity coefficient is taken into account, while refs. [1,22] discuss the variability of both coefficients. This form of the DPLE is also considered in the presented article.

It should be noted here that DPLE is also applicable in the field of bioheat transfer, e.g., [23,24]. Due to its specific structure, the biological tissue is characterized by significant delay times (from a few to several seconds). Therefore, equations that include lag times better approximate the course of thermal processes in the tissue domain than do the commonly known classical Pennes equation. The data reported in the literature regarding the value of delay times in biological tissues are very diverse, which is probably due to the consideration of different types of tissues and the difficulties associated with this type of experimental research [25].

As a rule, the practical problems related to the modeling of micro-scale heat transfer based on models that include delay times are solved using the numerical methods, but effective analytical or semi-analytical solutions can be also found (see: [26,27,28,29,30,31,32]). Another article [26] presents semi-analytical solutions for two-layer domains oriented in Cartesian, cylindrical and spherical coordinate systems (1D task), while a boundary condition with thermal resistance was assumed at the interface of the layers. In [27], the analytical solution to a dual-phase lag equation under the assumption that τ_T > τ_q using the separation of the variables and the Green’s function methods was obtained and presented. In [28], the dual-phase lag equation for the 3D problem is solved using the Adomian decomposition method, while ref. [29] presents the analytical solution to the 2D problem using the Green’s function method. In turn, in [30], a finite integral transform was used to solve the DPL equation describing bioheat transfer in heated tissues. Another paper [31] presents an analysis of a cracked strip irradiated by an ultrafast laser beam, while in [32], an analytical solution is obtained using the Duhamel’s theorem and the finite integral transform.

These are usually 1D tasks supplemented with uncomplicated boundary and initial conditions.

Numerical solutions to the equations with lag times (in particular for the tasks relating to more practical problems) are most often solved using the different variants of the finite difference method (FDM) or the control volume method (in some cases, the systems of equations that need to be solved in the successive time steps are identical). At the stage of numerical modeling, explicit, implicit and mixed schemes are applied, e.g., [1,11,33,34]. In particular, in [1], the implicit scheme of FDM supplemented by Gauss–Seidl iterative procedure has been used. In [11], the compact alternating direction method is used, while in [33] the FDM solution for double-layered thin film is presented. In turn, paper [34] presents the corrective smoothed particle method.

In the literature, one can find the numerical solutions to DPLE based on the finite element method [35,36] or the boundary element method [37].

The sensitivity analysis discussed in the paper presented concerns the changes of non-steady temperature field in the metal film domain due to perturbations in lag times. The adopted values of delay times are of fundamental importance to the quality of the numerical solution to the problem. If disturbances in delay times cause small changes in the calculation results, the problem becomes insignificant. Otherwise, the delay times should be carefully determined (measured). The remaining parameters occurring in DPLE are much easier to measure (determine); therefore, these two key delay times related to DPL were subjected to sensitivity analysis. Most often, the influence of parameters on temperatures is analyzed by taking their different values and showing temperature distributions. Sensitivity analysis is a more universal tool—it allows us to estimate the influence of parameter changes on temperature changes without the need to perform a very large number of computations. It should be noted that sensitivity analysis is usually performed under the assumption that all thermophysical parameters appearing in the equation are constant. In this paper, we assume that the thermal conductivity and the volumetric specific heat are temperature-dependent, which makes the task more difficult.

As is well known, the definition of the sensitivity function U_k (X, t) is as follows:

$$U_k(X,t)=\underset{\Delta p_k\to 0}{\lim }\frac{T(X,t,p_1,\dots ,p_k+\Delta p_k,\dots ,p_n)-T(X,t,p_1,\dots ,p_k,\dots ,p_n)}{\Delta p_k}=\frac{\partial T(X,t)}{\partial p_k}$$

(2)

Thus, the sensitivity model can be constructed using the differentiation of DPLE and the appropriate boundary-initial conditions with respect to the parameter p_k. This rather obvious approach is called a direct method—e.g., [38,39,40].

To avoid the need to create a mathematical model of sensitivity analysis, in engineering practice, Equation (2) is often used directly, replacing the actual sensitivity value with a difference quotient approximating the first derivative. Thus, the basic problem is solved twice, introducing small disturbances to the parameter in question, and it allows us to approximate the local and temporary sensitivity values at the selected points from the domain considered. Of course, this approach is a black box and does not allow for more detailed research related to sensitivity analysis.

2. Materials and Methods

2.1. Mathematical Description of the Process

The basic version of DPLE results from the expansion of the both sides of Equation (1) into the Taylor series with the accuracy to the first derivative, as follows:

$$q\left(X,t\right)+{\tau }_q\frac{\partial q\left(X,t\right)}{\partial t}=-\lambda (T)\mathrm{grad}T\left(X,t\right)-{\tau }_T\lambda (T)\frac{\partial \left[\mathrm{grad}T\left(X,t\right)\right]}{\partial t}$$

(3)

Next, the well-known Fourier equation is used, as follows:

$$\begin{array}{c}C\left(T\right)\frac{\partial T\left(X,t\right)}{\partial t}=-\mathrm{div}q\left(X,t\right)+Q\left(X,t\right)\end{array}$$

(4)

where C(T) is a volumetric specific heat and Q(X, t) is a capacity of internal heat sources.

After appropriate mathematical transformations and taking into account the variability of thermophysical parameters with temperature (detailed considerations can be found in [1]), the following form of an equation with two delay times is obtained:

$$\begin{array}{c}\begin{array}{l}C(T)\left(\frac{\partial T(X,t)}{\partial t}+{\tau }_q\frac{{\partial }^{2}T(X,t)}{\partial t^2}\right)+{\tau }_q\frac{\mathrm{d}C(T)}{\mathrm{d}T}{\left(\frac{\partial T(X,t)}{\partial t}\right)}^2=\\ \mathrm{div}\left[\lambda (T)\mathrm{grad}T(X,t)\right]+{\tau }_T\mathrm{div}\left[\lambda (T)\mathrm{grad}\left(\frac{\partial T(X,t)}{\partial t}\right)\right]+Q(X,t)+{\tau }_q\frac{\partial Q(X,t)}{\partial t}\end{array}\end{array}$$

(5)

In the case under consideration, the efficiency of internal sources Q(X, t) corresponds to the heat source resulting from the laser interaction. This approach is very often used in the problems presented here and is called volumetric source formulation [41].

The geometry of the laser beam and the resulting form of the function describing the efficiency of internal heat sources suggest treating the task as axially symmetrical and orient the domain in a cylindrical coordinate system (Figure 1). Thus, a certain cylindrical sub-domain is cut out from the metal layer. The dimensions of the domain are large enough to assume the adiabatic conditions on the external boundary.

The differential operator containing the geometrical coordinates for an axisymmetric problem is of the form

$$\begin{array}{l}\begin{array}{c}\mathrm{div}\left[\lambda (T)\mathrm{grad}T\right]+\end{array}{\tau }_T\mathrm{div}\left[\lambda (T)\mathrm{grad}\left(\frac{\partial T}{\partial t}\right)\right]=\\ \frac{1}{r}\frac{\partial }{\partial r}\left[r\lambda (T)\frac{\partial T}{\partial r}\right]+\frac{\partial }{\partial z}\left[\lambda (T)\frac{\partial T}{\partial z}\right]+\\ {\tau }_T\left\{\frac{1}{r}\frac{\partial }{\partial r}\left[r\lambda (T)\frac{\partial }{\partial r}\left(\frac{\partial T}{\partial t}\right)\right]+\frac{\partial }{\partial z}\left[\lambda (T)\frac{\partial }{\partial z}\left(\frac{\partial T}{\partial t}\right)\right]\right\}=\\ \frac{1}{r}\lambda (T)\left[\frac{\partial T}{\partial r}+{\tau }_T\frac{\partial }{\partial r}\left(\frac{\partial T}{\partial t}\right)\right]+\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}\frac{\partial T}{\partial r}\left[\frac{\partial T}{\partial r}+{\tau }_T\frac{\partial }{\partial r}\left(\frac{\partial T}{\partial t}\right)\right]+\\ \lambda (T)\left[\frac{{\partial }^{2}T}{\partial r^2}+{\tau }_T\frac{{\partial }^2}{\partial r^2}\left(\frac{\partial T}{\partial t}\right)\right]+\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}\frac{\partial T}{\partial z}\left[\frac{\partial T}{\partial z}+{\tau }_T\frac{\partial }{\partial z}\left(\frac{\partial T}{\partial t}\right)\right]+\\ \lambda (T)\left[\frac{{\partial }^{2}T}{\partial z^2}+{\tau }_T\frac{{\partial }^2}{\partial z^2}\left(\frac{\partial T}{\partial t}\right)\right]\end{array}$$

(6)

where r and z are the radial and axial co-ordinates.

The adiabatic Neumann boundary condition in the case of DPLE (see: [22]) is of the form

$$-\lambda (T)\left(\mathsf{n}\cdot \mathrm{grad}T(r,z,t)+{\tau }_T\mathsf{n}\cdot \mathrm{grad}\frac{\partial T(r,z,t)}{\partial t}\right)=0$$

(7)

where n is the normal outward vector. In the case under consideration, the normal derivative at the points on the side surface of the cylinder corresponds to derivative $\frac{\partial T(r,z,t)}{\partial r}$, and at the upper and lower bases, it corresponds to derivative $\pm \frac{\partial T(r,z,t)}{\partial z}$.

In this paper, the mathematical formula determining the intensity of heat source Q(r, z, t) is accepted in the following form [33,42]:

$$Q(r,z,t)=\sqrt{\frac{4\ln 2}{\pi }}(1-R)\frac{I_0}{\delta t_p}\exp \left[-\frac{r^2}{r_D^2}-\frac{z}{\delta }-4\ln 2\frac{{(t-2t_p)}^2}{t_p^2}\right]$$

(8)

where I₀ [J/m²] is the laser intensity, t_p [s] is the characteristic time of laser pulse, δ [m] is the optical penetration depth, R is the reflectivity of the irradiated surface, and r_D [m] is the laser beam radius.

It should be noted that the way the laser effect is taken into account depends on the material under consideration. For highly absorbed materials, laser heating is described by the Neumann boundary condition, and in Equation (5), the internal source function Q is equal to zero.

For strongly scattering materials (such as gold and nickel considered here), laser heating is included in the source component Q, and a thermal-insulation condition is assumed on the heated surface.

In turn, in [41], the other form of the boundary condition on the laser-heated surface is given. It is quite complicated in form from the mathematical point of view, and the adequate formula is similar to the Robin boundary condition. The condition discussed in [41] contains a certain parameter varying from 0 to 1, the numerical value of which is not easy to determine.

The initial conditions are also known, as follows:

$$t=0:T(r,z,t)=T_p,\frac{\partial T(r,z,t)}{\partial t}=0$$

(9)

where T_p is the initial temperature of the domain.

2.2. Sensitivity Models

As mentioned in the introduction, a direct method was used to create sensitivity-analysis equations and involves differentiating the basic equation and boundary-initial conditions with respect to the parameter under consideration [38,39].

Below, subsequent transformations leading to the final form of the sensitivity equations with respect to delay times are presented.

Thus, the energy equation must be differentiated on both sides with respect to the relaxation time, as follows:

$$\begin{array}{l}\begin{array}{c}\begin{array}{l}\frac{\partial C(T)}{\partial {\tau }_q}\frac{\partial T}{\partial t}+C(T)\frac{\partial }{\partial {\tau }_q}\left(\frac{\partial T}{\partial t}\right)+C(T)\frac{{\partial }^{2}T}{\partial t^2}+{\tau }_q\frac{\partial C(T)}{\partial {\tau }_q}\frac{{\partial }^{2}T}{\partial t^2}+{\tau }_{q}C(T)\frac{\partial }{\partial {\tau }_q}\left(\frac{{\partial }^{2}T}{\partial t^2}\right)+\\ \frac{\mathrm{d}C(T)}{\mathrm{d}T}{\left(\frac{\partial T}{\partial t}\right)}^2+{\tau }_q\frac{\partial }{\partial {\tau }_q}\left(\frac{\mathrm{d}C(T)}{\mathrm{d}T}\right){\left(\frac{\partial T}{\partial t}\right)}^2+{\tau }_q\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial }{\partial {\tau }_q}\left({\left(\frac{\partial T}{\partial t}\right)}^2\right)=\end{array}\end{array}\\ \mathrm{div}\left[\frac{\partial \lambda (T)}{\partial {\tau }_q}\mathrm{grad}T\right]+\mathrm{div}\left[\lambda (T)\mathrm{grad}\frac{\partial T}{\partial {\tau }_q}\right]+{\tau }_T\mathrm{div}\left[\frac{\partial \lambda (T)}{\partial {\tau }_q}\mathrm{grad}\left(\frac{\partial T}{\partial t}\right)\right]+\\ {\tau }_T\mathrm{div}\left[\lambda (T)\mathrm{grad}\frac{\partial }{\partial {\tau }_q}\left(\frac{\partial T}{\partial t}\right)\right]+\frac{\partial Q}{\partial t}\end{array}$$

(10)

Because

$$\begin{array}{c}\frac{\partial C(T)}{\partial {\tau }_q}=\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial T}{\partial {\tau }_q}\end{array}$$

(11)

consequently,

$$\begin{array}{c}\begin{array}{l}\frac{\partial }{\partial {\tau }_q}\left(\frac{\mathrm{d}C(T)}{\mathrm{d}T}\right)=\frac{\mathrm{d}}{\mathrm{d}T}\left(\frac{\partial C(T)}{\partial {\tau }_q}\right)=\frac{\mathrm{d}}{\mathrm{d}T}\left(\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial T}{\partial {\tau }_q}\right)=\frac{{\mathrm{d}}^{2}C(T)}{\mathrm{d}{T}^2}\frac{\partial T}{\partial {\tau }_q}+\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\mathrm{d}}{\mathrm{d}T}\left(\frac{\partial T}{\partial {\tau }_q}\right)=\\ \frac{{\mathrm{d}}^{2}C(T)}{\mathrm{d}{T}^2}\frac{\partial T}{\partial {\tau }_q}+\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial }{\partial {\tau }_q}\left(\frac{\mathrm{d}T}{\mathrm{d}T}\right)=\frac{{\mathrm{d}}^{2}C(T)}{\mathrm{d}{T}^2}\frac{\partial T}{\partial {\tau }_q}+\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial }{\partial {\tau }_q}\left(\frac{\mathrm{d}T}{\mathrm{d}T}\right)=\frac{{\mathrm{d}}^{2}C(T)}{\mathrm{d}{T}^2}\frac{\partial T}{\partial {\tau }_q}\end{array}\end{array}$$

(12)

and

$$\begin{array}{c}\frac{\partial }{\partial {\tau }_q}\left({\left(\frac{\partial T}{\partial t}\right)}^2\right)=\frac{\partial }{\partial {\tau }_q}\left(\frac{\partial T}{\partial t}\frac{\partial T}{\partial t}\right)=\frac{\partial }{\partial t}\left(\frac{\partial T}{\partial {\tau }_q}\right)\frac{\partial T}{\partial t}+\frac{\partial T}{\partial t}\frac{\partial }{\partial t}\left(\frac{\partial T}{\partial {\tau }_q}\right)=2\frac{\partial }{\partial t}\left(\frac{\partial T}{\partial {\tau }_q}\right)\frac{\partial T}{\partial t}\end{array}$$

(13)

Similarly,

$$\begin{array}{c}\frac{\partial \lambda (T)}{\partial {\tau }_q}=\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}\frac{\partial T}{\partial {\tau }_q}\end{array}$$

(14)

that is,

$$\begin{array}{l}\begin{array}{c}\begin{array}{l}\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial T}{\partial {\tau }_q}\frac{\partial T}{\partial t}+C(T)\frac{\partial }{\partial t}\left(\frac{\partial T}{\partial {\tau }_q}\right)+C(T)\frac{{\partial }^{2}T}{\partial t^2}+{\tau }_q\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial T}{\partial {\tau }_q}\frac{{\partial }^{2}T}{\partial t^2}+{\tau }_{q}C(T)\frac{{\partial }^2}{\partial t^2}\left(\frac{\partial T}{\partial {\tau }_q}\right)+\\ \frac{\mathrm{d}C(T)}{\mathrm{d}T}{\left(\frac{\partial T}{\partial t}\right)}^2+{\tau }_q\frac{{\mathrm{d}}^{2}C(T)}{\mathrm{d}{T}^2}\frac{\partial T}{\partial {\tau }_q}{\left(\frac{\partial T}{\partial t}\right)}^2+2{\tau }_q\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial }{\partial t}\left(\frac{\partial T}{\partial {\tau }_q}\right)\frac{\partial T}{\partial t}=\end{array}\end{array}\\ \mathrm{div}\left[\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}\frac{\partial T}{\partial {\tau }_q}\mathrm{grad}T\right]+\mathrm{div}\left[\lambda (T)\mathrm{grad}\frac{\partial T}{\partial {\tau }_q}\right]+{\tau }_T\mathrm{div}\left[\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}\frac{\partial T}{\partial {\tau }_q}\mathrm{grad}\left(\frac{\partial T}{\partial t}\right)\right]+\\ {\tau }_T\mathrm{div}\left[\lambda (T)\mathrm{grad}\frac{\partial }{\partial t}\left(\frac{\partial T}{\partial {\tau }_q}\right)\right]+\frac{\partial Q}{\partial t}\end{array}$$

(15)

or

$$\begin{array}{l}C(T)\left(\frac{\partial U}{\partial t}+{\tau }_q\frac{{\partial }^{2}U}{\partial t^2}\right)+\frac{\mathrm{d}C(T)}{\mathrm{d}T}U\left(\frac{\partial T}{\partial t}+{\tau }_q\frac{{\partial }^{2}T}{\partial t^2}\right)+\\ +C(T)\frac{{\partial }^{2}T}{\partial t^2}+\left(\frac{\mathrm{d}C(T)}{\mathrm{d}T}+{\tau }_q\frac{{\mathrm{d}}^{2}C(T)}{\mathrm{d}{T}^2}U\right){\left(\frac{\partial T}{\partial t}\right)}^2+2{\tau }_q\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial U}{\partial t}\frac{\partial T}{\partial t}=\\ \mathrm{div}\left[\lambda (T)\mathrm{grad}U\right]+{\tau }_T\mathrm{div}\left[\lambda (T)\mathrm{grad}\frac{\partial U}{\partial t}\right]+\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}\mathrm{div}\left(U\mathrm{grad}T\right)+\\ {\tau }_T\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}\mathrm{div}\left[U\mathrm{grad}\left(\frac{\partial T}{\partial t}\right)\right]+\frac{\partial Q}{\partial t}\end{array}$$

(16)

where U = ∂T/∂τ_q is the sensitivity with respect to relaxation time.

The adiabatic Neumann boundary condition (7) takes the form

$$-\frac{\partial \lambda (T)}{\partial {\tau }_q}\left(\mathsf{n}\cdot \mathrm{grad}T+{\tau }_T\mathsf{n}\cdot \mathrm{grad}\frac{\partial T}{\partial t}\right)-\lambda (T)\left(\mathsf{n}\cdot \mathrm{grad}\frac{\partial T}{\partial {\tau }_q}+{\tau }_T\mathsf{n}\cdot \mathrm{grad}\frac{\partial }{\partial {\tau }_q}\left(\frac{\partial T}{\partial t}\right)\right)=0$$

(17)

or

$$-\lambda (T)\left(\mathsf{n}\cdot \mathrm{grad}U+{\tau }_T\mathsf{n}\cdot \mathrm{grad}\frac{\partial U}{\partial t}\right)=\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}U\left(\mathsf{n}\cdot \mathrm{grad}T+{\tau }_T\mathsf{n}\cdot \mathrm{grad}\frac{\partial T}{\partial t}\right)$$

(18)

The initial condition of the sensitivity model is as follows:

$$t=0:\frac{\partial T}{\partial {\tau }_q}=\frac{\partial T_p}{\partial {\tau }_q},\frac{\partial }{\partial {\tau }_q}\left(\frac{\partial T}{\partial t}\right)=0$$

(19)

this means

$$t=0:U=0,\frac{\partial U}{\partial t}=0$$

(20)

A similar procedure is followed when creating the sensitivity equation with respect to the thermalization time. Thus, we differentiate the energy equation with respect to τ_T, as follows:

$$\begin{array}{c}\begin{array}{l}\frac{\partial C(T)}{\partial {\tau }_T}\frac{\partial T}{\partial t}+C(T)\frac{\partial }{\partial {\tau }_T}\left(\frac{\partial T}{\partial t}\right)+{\tau }_q\frac{\partial C(T)}{\partial {\tau }_T}\frac{{\partial }^{2}T}{\partial t^2}+{\tau }_{q}C(T)\frac{\partial }{\partial {\tau }_T}\left(\frac{{\partial }^{2}T}{\partial t^2}\right)+\\ {\tau }_q\frac{\partial }{\partial {\tau }_T}\left(\frac{\mathrm{d}C(T)}{\mathrm{d}T}\right){\left(\frac{\partial T}{\partial t}\right)}^2+{\tau }_q\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial }{\partial {\tau }_T}\left({\left(\frac{\partial T}{\partial t}\right)}^2\right)=\\ \mathrm{div}\left[\frac{\partial \lambda (T)}{\partial {\tau }_T}\mathrm{grad}T\right]+\mathrm{div}\left[\lambda (T)\mathrm{grad}\left(\frac{\partial T}{\partial {\tau }_T}\right)\right]+\mathrm{div}\left[\lambda (T)\mathrm{grad}\left(\frac{\partial T}{\partial t}\right)\right]+\\ {\tau }_T\mathrm{div}\left[\frac{\partial \lambda (T)}{\partial {\tau }_T}\mathrm{grad}\left(\frac{\partial T}{\partial t}\right)\right]+{\tau }_T\mathrm{div}\left[\lambda (T)\mathrm{grad}\frac{\partial }{\partial {\tau }_T}\left(\frac{\partial T}{\partial t}\right)\right]\end{array}\end{array}$$

(21)

or

$$\begin{array}{c}\begin{array}{l}\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial T}{\partial {\tau }_T}\frac{\partial T}{\partial t}+C(T)\frac{\partial }{\partial t}\left(\frac{\partial T}{\partial {\tau }_T}\right)+{\tau }_q\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial T}{\partial {\tau }_T}\frac{{\partial }^{2}T}{\partial t^2}+{\tau }_{q}C(T)\frac{{\partial }^2}{\partial t^2}\left(\frac{\partial T}{\partial {\tau }_T}\right)+\\ {\tau }_q\frac{{\mathrm{d}}^{2}C(T)}{\mathrm{d}{T}^2}\frac{\partial T}{\partial {\tau }_T}{\left(\frac{\partial T}{\partial t}\right)}^2+2{\tau }_q\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial }{\partial t}\left(\frac{\partial T}{\partial {\tau }_T}\right)\frac{\partial T}{\partial t}=\\ \mathrm{div}\left[\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}\frac{\partial T}{\partial {\tau }_T}\mathrm{grad}T\right]+\mathrm{div}\left[\lambda (T)\mathrm{grad}\left(\frac{\partial T}{\partial {\tau }_T}\right)\right]+\mathrm{div}\left[\lambda (T)\mathrm{grad}\left(\frac{\partial T}{\partial t}\right)\right]+\\ {\tau }_T\mathrm{div}\left[\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}\frac{\partial T}{\partial {\tau }_T}\mathrm{grad}\left(\frac{\partial T}{\partial t}\right)\right]+{\tau }_T\mathrm{div}\left[\lambda (T)\mathrm{grad}\frac{\partial }{\partial t}\left(\frac{\partial T}{\partial {\tau }_T}\right)\right]\end{array}\end{array}$$

(22)

and finally

$$\begin{array}{c}\begin{array}{l}C(T)\left(\frac{\partial W}{\partial t}+{\tau }_q\frac{{\partial }^{2}W}{\partial t^2}\right)+\frac{\mathrm{d}C(T)}{\mathrm{d}T}W\left(\frac{\partial T}{\partial t}+{\tau }_q\frac{{\partial }^{2}T}{\partial t^2}\right)+\\ {\tau }_q\frac{{\mathrm{d}}^{2}C(T)}{\mathrm{d}{T}^2}W{\left(\frac{\partial T}{\partial t}\right)}^2+2{\tau }_q\frac{\mathrm{d}C(T)}{\mathrm{d}T}\frac{\partial W}{\partial t}\frac{\partial T}{\partial t}=\\ \mathrm{div}\left[\lambda (T)\mathrm{grad}W\right]+{\tau }_T\mathrm{div}\left[\lambda (T)\mathrm{grad}\frac{\partial W}{\partial t}\right]+\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}\mathrm{div}\left(W\mathrm{grad}T\right)\\ +\mathrm{div}\left[\lambda (T)\mathrm{grad}\left(\frac{\partial T}{\partial t}\right)\right]+{\tau }_T\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}\mathrm{div}\left[W\mathrm{grad}\left(\frac{\partial T}{\partial t}\right)\right]\end{array}\end{array}$$

(23)

where W = ∂T/∂τ_T is the sensitivity with respect to thermalization time.

The adiabatic Neumann boundary condition complementary to Equation (23)

$$\begin{array}{l}-\frac{\partial \lambda (T)}{\partial {\tau }_T}\left(\mathsf{n}\cdot \mathrm{grad}T+{\tau }_T\mathsf{n}\cdot \mathrm{grad}\frac{\partial T}{\partial t}\right)-\\ \lambda (T)\left(\mathsf{n}\cdot \mathrm{grad}\frac{\partial T}{\partial {\tau }_T}+\mathsf{n}\cdot \mathrm{grad}\frac{\partial T}{\partial t}+{\tau }_T\mathsf{n}\cdot \mathrm{grad}\frac{\partial }{\partial {\tau }_T}\left(\frac{\partial T}{\partial t}\right)\right)=0\end{array}$$

(24)

can be written in the form

$$\begin{array}{l}-\lambda (T)\left(\mathsf{n}\cdot \mathrm{grad}W+{\tau }_T\mathsf{n}\cdot \mathrm{grad}\frac{\partial W}{\partial t}\right)=\\ \frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}W\left(\mathsf{n}\cdot \mathrm{grad}T+{\tau }_T\mathsf{n}\cdot \mathrm{grad}\frac{\partial T}{\partial t}\right)+\lambda (T)\mathsf{n}\cdot \mathrm{grad}\frac{\partial T}{\partial t}\end{array}$$

(25)

while the initial condition is the following:

$$t=0:W=0,\frac{\partial W}{\partial t}=0$$

(26)

2.3. Method of Solution

As mentioned previously, the numerical algorithm based on the implicit scheme of the finite-difference method has been used. The domain is covered by the regular mesh with the step $\Delta r=\Delta z=h$. Inside the domain the five-points stars have been distinguished, namely P_i,j, P_i+1,j, P_i−1,j, P_i,j+1, P_i,j−1,j while i = 1, …, m−1, j = 1, …, n − 1. The boundary nodes corresponds to i = 0, i = m, j = 0 and j = n. Time step $\Delta t$ is assumed to be constant. The succesive three time levels are denoted as f − 2, f − 1 and f. For the above notations, the differential form of DPLE is of the following form (detailed considerations can be found in [1]):

$$\begin{array}{l}{C}_{i,j}^{f-1}\frac{T_{i,j}^f-T_{i,j}^{f-1}}{\Delta t}+{\tau }_q{\left(\frac{\mathrm{d}C\left(T\right)}{\mathrm{d}T}\right)}_{i,j}^{f-1}{\left(\frac{T_{i,j}^{f-1}-T_{i,j}^{f-2}}{\Delta t}\right)}^2+C_{i,j}^{f-1}{\tau }_q\frac{T_{i,j}^f-2T_{i,j}^{f-1}+T_{i,j}^{f-2}}{{\left(\Delta t\right)}^2}=\\ \frac{{\lambda }_{i,j}^{f-1}\left(\Delta t+{\tau }_T\right)}{h^2\Delta t}\left(1-\frac{h}{2r_{i,j}}\right)T_{i,j-1}^f-4\frac{{\lambda }_{i,j}^{f-1}\left(\Delta t+{\tau }_T\right)}{h^2\Delta t}{T}_{i,j}^f+\\ \frac{{\lambda }_{i,j}^{f-1}\left(\Delta t+{\tau }_T\right)}{h^2\Delta t}\left(1+\frac{h}{2r_{i,j}}\right)T_{i,j+1}^f+\frac{{\lambda }_{i,j}^{f-1}\left(\Delta t+{\tau }_T\right)}{h^2\Delta t}\left(T_{i+1,j}^f+T_{i-1,j}^f\right)+R_{i,j}^{f-1}+\\ Q_{i,j}^f+{\tau }_q{\left(\frac{\partial Q}{\partial t}\right)}_{i,j}^f\end{array}$$

(27)

where

$$\begin{array}{l}{R}_{i,j}^{f-1}=-\frac{{\lambda }_{i,j}^{f-1}{\tau }_T}{h^2\Delta t}\left(T_{i,j-1}^{f-1}+T_{i,j+1}^{f-1}+T_{i-1,j}^{f-1}+T_{i+1,j}^{f-1}-4T_{i,j}^{f-1}\right)-\frac{{\lambda }_{i,j}^{f-1}{\tau }_T}{2h{r}_{i,j}\Delta t}\left(T_{i,j+1}^{f-1}-T_{i,j-1}^{f-1}\right)+\\ \frac{1}{4h^2\Delta t}{\left(\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}\right)}_{i,j}^{f-1}\left(T_{i,j+1}^{f-1}-T_{i,j-1}^{f-1}\right)\left[\left(\Delta t+{\tau }_T\right)\left(T_{i,j+1}^{f-1}-T_{i,j-1}^{f-1}\right)-{\tau }_T\left(T_{i,j+1}^{f-2}-T_{i,j-1}^{f-2}\right)\right]+\\ \frac{1}{4h^2\Delta t}{\left(\frac{\mathrm{d}\lambda (T)}{\mathrm{d}T}\right)}_{i,j}^{f-1}\left(T_{i+1,j}^{f-1}-T_{i-1,j}^{f-1}\right)\left[\left(\Delta t+{\tau }_T\right)\left(T_{i+1,j}^{f-1}-T_{i-1,j}^{f-1}\right)-{\tau }_T\left(T_{i+1,j}^{f-2}-T_{i-1,j}^{f-2}\right)\right]\end{array}$$

(28)

After not-too-complex mathematical manipulations, one obtains the following final form of Equation (27):

$$\begin{array}{l}{T}_{i,j}^f=\frac{{\lambda }_{i,j}^{f-1}\left(\Delta t+{\tau }_T\right)}{A_{i,j}^{f-1}{h}^2\Delta t}\left(1-\frac{h}{2r_{i,j}}\right)T_{i,j-1}^f+\frac{{\lambda }_{i,j}^{f-1}\left(\Delta t+{\tau }_T\right)}{A_{i,j}^{f-1}{h}^2\Delta t}\left(1+\frac{h}{2r_{i,j}}\right)T_{i,j+1}^f+\\ \frac{{\lambda }_{i,j}^{f-1}\left(\Delta t+{\tau }_T\right)}{A_{i,j}^{f-1}{h}^2\Delta t}\left(T_{i+1,j}^f+T_{i-1,j}^f\right)+\frac{B_{i,j}^{f-1}}{A_{i,j}^{f-1}}\end{array}$$

(29)

where

$$\begin{array}{l}{A}_{i,j}^{f-1}=\frac{C_{i,j}^{f-1}\left(\Delta t+{\tau }_q\right)}{{\left(\Delta t\right)}^2}+4\frac{{\lambda }_{i,j}^{f-1}\left(\Delta t+{\tau }_T\right)}{h^2\Delta t}\\ B_{i,j}^{f-1}=R_{i,j}^{f-1}+\frac{C_{i,j}^{f-1}\left(\Delta t+2{\tau }_q\right)}{{\left(\Delta t\right)}^2}{T}_{i,j}^{f-1}-\frac{C_{i,j}^{f-1}{\tau }_q}{{\left(\Delta t\right)}^2}{T}_{i,j}^{f-2}-\\ {\tau }_q{\left(\frac{\mathrm{d}C\left(T\right)}{\mathrm{d}T}\right)}_{i,j}^{f-1}{\left(\frac{T_{i,j}^{f-1}-T_{i,j}^{f-2}}{\Delta t}\right)}^2+Q_{i,j}^f+{\tau }_q{\left(\frac{\partial Q}{\partial t}\right)}_{i,j}^f\end{array}$$

(30)

The numerical approximation of the boundary conditions is quite simple and will not be presented here [1].

The differential approximation of subsequent components of the energy equation was adopted in such a way that the final equations of the finite difference method algorithm are the linear ones. This system of equations was solved using the Gauss–Seidel iterative procedure.

As can be seen, the sensitivity equations with respect to delay times have a form quite similar to that of the basic energy equation. Thus, the numerous fragments of the program used for the DPLE numerical solution can be applied to construct the numerical algorithm of the solution to sensitivity equations, and the final form of the finite difference method equations does not differ significantly from Equation (29). Of course, both parts of the program must be solved step-by-step simultaneously because determining the instantaneous sensitivity field requires knowledge of the temperature field at the same moment in time.

3. Results of Computations and Their Discussion

We consider the cylindrical domain R₀ = 100 nm, Z = 100 nm subjected to the laser pulse, with the initial temperature of domain equal to T_p = 300 K. The first numerical solutions correspond to the thin metal film made from gold, while the second set correspond to the film made from nickel. The functions describing the temperature-dependent thermal conductivity and volumetric specific heat of gold can be found in [43], and they are of the forms

$$\lambda (T)=320.973-0.0111T-2.747\cdot {10}^{-5}{T}^2-4.048\cdot {10}^{-9}{T}^3[\mathrm{W}/(\mathrm{mK})]$$

(31)

and

$$\begin{array}{l}C(T)=19300(105.1+0.2941T-8.731\cdot {10}^{-4}{T}^2+1.787\cdot {10}^{-6}{T}^3-\\ 7.051\cdot {10}^{-10}{T}^4+1.538\cdot {10}^{-13}{T}^5)[{\mathrm{J}/(\mathrm{m}}^3\mathrm{K})]\end{array}$$

(32)

while the parameters of nickel [44] are described by

$$\lambda (T)[\mathrm{W}/(\mathrm{mK})]=\left\{\begin{array}{l}114.43-0.0082T,T\le 600\mathrm{K}?\\ 50.47+0.021T,600\mathrm{K}<T\le 1726\mathrm{K}\end{array}\right.$$

(33)

and

$$C(T)[{\mathrm{J}/(\mathrm{m}}^3\mathrm{K})]=8900\left\{\begin{array}{l}-295.95+4.95T-0.0096T^2+9.46\cdot {10}^{-6}{T}^3-\\ 4.36\cdot {10}^{-9}{T}^4+7.6\cdot {10}^{-13}{T}^5,T\le 1400\mathrm{K}?\\ 616.56,1400\mathrm{K}<T\le 1726\mathrm{K}?\end{array}\right.$$

(34)

These relationships were used for temperatures lower than the melting point (gold, 1336 K; nickel, 1726 K).

The relaxation and thermalization times for gold are equal to 8.5 ps and 90 ps [43], while for nickel, they are 0.82 ps and 10 ps, respectively [44].

In the cases presented below, half of the cylinder cross-section was covered with the regular geometric mesh containing $100\times 100$ nodes (mesh step h = 1 nm), The analysis of numerical computation results confirmed that such a mesh density is fully satisfactory (meshes $50\times 50$ and $200\times 200$ have been also tested). A time step of Δt = 0.0005 ps has been assumed. The relatively small size of the time step was intended to provide a good approximation of the temporary efficiency of an internal heat source Q(r, z, t).

The example of temperature distribution in the domain considered is shown in Figure 2. The presented results confirm, among other points, the correctness of the adiabatic boundary conditions assumption on the external surfaces of the domain.

Figure 3 illustrates the temperature history at points A(0, 1 nm), B(0, 20 nm) and C(0, 40 nm) from the interior of domain considered. The heating/cooling curves do not differ qualitatively from typical ones for intensive laser heating of thin metal films.

It is assumed below that the disturbances in delay times correspond to 5% of the base value, i.e., $\Delta {\tau }_q=0.05{\tau }_q,\Delta {\tau }_T=0.05{\tau }_T$. The perturbations of the temporary temperature values at the selected points corresponding to these values are equal to

$${\mathrm{D}}_{1}T=U\Delta \tau {}_q,{\mathrm{D}}_{2}T=W\Delta \tau {}_T,\mathrm{D}T=\sqrt{{\left(U\Delta {\tau }_q\right)}^2+{\left(W\Delta {\tau }_T\right)}^2}$$

(35)

From a qualitative point of view, the curves illustrating temperature-field disturbances (Figure 4) are similar to the basic solution. However, a shift of the maximum values to the right can be observed. This shift is due to changes in the values of delay times (here their increases are positive). Disturbances resulting from a change in the thermalization time (Figure 5) have a similar course to those related to the relaxation time but are negative. In Figure 6, the temperature changes caused by a simultaneous perturbation of parameters τ_q and τ_T are shown.

It can be expected that reducing the efficiency of the internal heat source (i.e., reducing the power of the laser beam or increasing the time t_p) will result in significantly smaller disturbances in the temperature field. This rather obvious fact is confirmed, among others, by the results of one of the numerical simulations for which I₀ = 10,000 W/m² and t_p = 1 ps, as shown in Figure 7, Figure 8, Figure 9 and Figure 10.

It can be noted here that the changes in the temperature field (especially in the axis of the cylinder near the heated outer surface) are significant and increase with the increase in the heating intensity (c.f. Figure 6 and Figure 10).

The second series of numerical simulations concerns the laser heating of an ultra-thin nickel layer. Therefore, a material with much greater resistance to external thermal impacts has been considered. The geometrical parameters of the layer sub-domain are the same as in the previous examples. In Figure 11, the temperature distribution in the domain is shown, while Figure 12 illustrates the heating/cooling curves at points A and B. The curve at point C is not shown in Figure 12 because at this point, the temperature remains practically unchanged, while the sensitivity values are equal to zero.

Figure 13, Figure 14 and Figure 15 illustrate the changes in the temperature field found using the sensitivity analysis methods caused by 5% disturbances in the relaxation time, thermalization time and both delay times.

It can be seen that the changes in the temperature field caused by the relaxation time disturbance are slight. In fact, this time in the case of nickel is small (0.82 ps), which to some extent brings the model with delay times closer to the Fourier model. Additionally, taking into account the values of thermophysical parameters of this material, one can conclude that the heating/cooling rate is clearly lower than, for example, that found in the case of gold, and this also results in small changes to local and temporary temperature values.

To verify the results of the computations, the metal layer made of gold is considered (Z = 0.2 μm). The task is treated as a 1D problem. It results from the experimental data quoted in [45]). The boundary conditions on the external surfaces of the system are assumed to be adiabatic ones. The source function depends only on the coordinate z and

$$Q(z,t)=(1-R)\frac{I_0}{\delta t_p\sqrt{\pi }}\exp \left[-\frac{z}{\delta }-\frac{{(t-2t_p)}^2}{t_p^2}\right]$$

(36)

where R = 0.93, I₀ = 500 J/m², δ = 5 μm, t_p = 0.096 ps. The results of measurements presented in [45] are given in dimensionless form, namely [T(0) − T_p]/(T_max − T_p) at point 0, where T_max is the maximum temperature.

The application of the developed program to an axially symmetrical problem can be very easily adapted to a 1D task and was used in numerical calculations. In

Table 1. Comparison of the results.

Time [ps]	0.1	0.2	0.3	0.5	1.0	1.5	2.0	2.5	3.0
Experiment	0.11	0.40	1.0	0.65	0.38	0.24	0.20	0.12	0.10
Model	0..11	0.38	1.0	0.67	0.36	0.22	0.19	0.11	0.11

, the comparison of the numerical and experimental solutions is presented. It can be said that the comparison of the results is quite positive.

4. Conclusions

The topic of the work is an equation with two delay times that takes into account variable, temperature-dependent thermophysical parameters. This form of the equation is, of course, more complex than the classical DPLE. The research concerns the analysis of the sensitivity of the transient temperature field in the system in relation to delay times (relaxation and thermalization times). The theoretical part of the article gives the form of the basic energy equation and the appropriate boundary-initial conditions and presents the sensitivity equations with respect to the parameters considered. The sensitivity equations have been obtained by the direct method, i.e., by differentiating the basic equation and boundary-initial conditions with respect to the parameter under consideration. The second part of the work is devoted to the problems of the numerical solution t the formulated boundary-initial problems. In particular, an ultra-thin metal layer exposed to laser pulse is considered. A cylindrical sub-domain of material is cut out from the layer domain. Part of its upper surface is heated by a laser beam, and adiabatic conditions are assumed on the outer surfaces of the area. The laser effect has been taken into account with an additional component determining the efficiency of temporary and local internal heat sources (this approach is often used in the computations of this type). At the stage of numerical algorithm construction, the implicit scheme of the finite difference method (FDM) has been used. The differential operators approximating the derivatives in the basic equation and sensitivity equations were written in a way that ensures the linearity of the final system of equations, which is solved step-by-step using the Gauss–Seidel method. A detailed analysis of the results of numerical simulations and the conclusions resulting from them are presented in Section 3.

To the authors’ knowledge, the presented article is the first attempt to apply sensitivity-analysis methods to modeling thermal processes on the basis of a nonlinear equation with two delay times.

The following conclusions can be drawn from the results obtained:

-

taking into account the variability of thermophysical parameters (especially for higher temperatures) causes visible changes in the results of numerical simulations;

-

disturbances in the delay-time values clearly change the course of the heating/cooling process in the domain considered;

-

the sensitivity of the temperature field with respect to the delay times increases with the increase in laser intensity;

-

the sensitivity of the temperature field with respect to the delay times varies depending on the type of material and is greater when the metal has a higher mean conductivity coefficient and a lower mean volumetric specific heat.