Unconditional Stability of a Numerical Method for the Dual-Phase-Lag Equation

The stability properties of a numerical method for the dual-phase-lag (DPL) equation are analyzed. The DPL equation has been increasingly used to model microand nanoscale heat conduction in engineering and bioheat transfer problems. A discretization method for the DPL equation that could yield efficient numerical solutions of 3D problems has been previously proposed, but its stability properties were only suggested by numerical experiments. In this work, the amplificationmatrix of themethod is analyzed, and it is shown that its powers are uniformly bounded. As a result, the unconditional stability of the method is established.


Introduction
Non-Fourier heat conduction models have been increasingly used in recent years to model a variety of engineering and biological heat transfer problems (see, e.g., [1][2][3] and references therein).The dual-phase-lag model (DPL) proposed by Tzou [4,5] considers two phase-lags,   , corresponding to the heat flux vector q, and   , corresponding to the temperature gradient ∇, resulting in the non-Fourier heat conduction law q (r,  +   ) = −∇ (r,  +   ) , where  and r are the time and spatial coordinates and  is the thermal conductivity.
Considering first-order approximations in the phase-lags in (1), in the absence of heat sources and in one dimension (1D), the combination of (1) with the conservation of energy principle leads to the DPL equation where  is the thermal diffusivity and Δ is the Laplace operator.
As pointed out in [17], in the case of numerical finite difference methods conditions on the stability of the proposed methods are not always completely established.McDonough et al. [18] proposed a numerical solution procedure for the DPL equation that could provide the base for an efficient numerical solver for 3D problems.However, only a partial analysis of the stability of the method was provided.As indicated by the authors, a von Neumann analysis had been conducted, though not reported in their paper, but it only provided necessary, but not sufficient, criteria for stability.Also, the authors reported that no stability problems were observed in a wide set of numerical experiments, but, as they acknowledged, the unconditional stability of their method was only suggested, but not established.
The aim of this note was to prove the unconditional stability of the finite difference method proposed in [18], which will be accomplished by analyzing its amplification matrix (e.g., [19,20]).This type of von Neumann or Fourier stability analysis is directly valid for pure initial-value problems and for initial-boundary-value problems with periodic boundary conditions, and it can also be applied to general Dirichlet or Neumann boundary conditions by using appropriate periodic extensions, provided that consistent discretizations are used for the boundary conditions [21,Ch. 3].
In the next section we recall the method proposed in [18] and write the scheme in a form suitable for the posterior stability analysis, in the case where   ̸ = 0, letting the case   = 2 Mathematical Problems in Engineering 0 for a particular simplified analysis.Next, we will obtain the amplification matrix corresponding to the scheme, showing that their powers are uniformly bounded, thus establishing the unconditional stability of the method.Finally, the main conclusions are summarized.

Finite Difference Scheme
The numerical solution procedure for the DPL equation proposed by McDonough et al. [18] consists in applying first trapezoidal integration to (2) and then using the following finite difference approximations: where  = Δ > 0 and considering the spatial domain  ∈ [, ], with   uniformly spaced grid points, so that Δ = ℎ = ( − )/(  − 1) and  = 1, 2, . . .,   , and superscripts indicate time levels.The resulting finite difference scheme is globally firstorder accurate in time and second-order accurate in space.Assuming   ̸ = 0, the scheme can be written in the form where  = /ℎ 2 .The case   = 0 will be examined apart.

Unconditional Stability
2 stability of the three-level scheme (4) follows from the uniform boundedness in , for every 0 <  <  for a fixed , of the powers of the amplification matrix Ẽ() ([19, p. 69], [20, p. 65]): where For scheme (4), the coefficients   and   are as follows: Therefore, Introducing the parameter  = cos  − 1, (8) take the form Thus, writing one gets Finally, writing the amplification matrix Ẽ() can be written in the simple form We note the following bounds, which will be of use to obtain bounds on the eigenvalues of the amplification matrix.Since −1 ≤ cos  ≤ 1, it follows that −2 ≤  ≤ 0. Also, from (10), it follows that and, hence, Therefore, the following bounds for the parameters defined in (12) hold: Also, since  = (), there is a fixed  > 0 such that 1++ > 0 for every 0 <  < .
Next we compute the eigenvalues of the amplification matrix, as written in (13), whose characteristic equation is Thus, the eigenvalues of Ẽ() are given by writing  1 for + sign and  2 for − sign.In the case where there are two complex conjugate eigenvalues and also when there is a double root, from (17) one has that  1  2 =  < 1, and, therefore, When there are two real different eigenvalues, that is, when (1 +  + ) 2 − 4 > 0, one has  1 >  2 > 0, and also it follows that  1 ≤ 1, since  < 0 and Therefore, the following bound on the spectral radius of the amplification matrix holds: which shows that the scheme satisfies the von Neumann necessary condition for stability, as stated in [18].Moreover, it also holds that Ẽ() has at least one eigenvalue with modulus strictly lower than 1, which will be the complementary result used to prove the unconditional stability of the scheme.We will next follow similar arguments to those used in [20, pp. 66-68].First we note that the sum of the absolute values of the elements of the amplification matrix is bounded: Then, since Ẽ() is a (2 × 2) matrix, there is a unitary matrix  such that  =  * Ẽ() is the triangular matrix where Hence,  2 stability of scheme ( 4) is equivalent to the uniform boundedness of the powers of the triangular matrix : where Since max(| 1 |, | 2 |) ≤ 1, it is only necessary to prove the boundedness of   ( 1 ,  2 ).When the eigenvalues are complex conjugate or real double, there is  < 1 such that Also, when there are two different real eigenvalues, it holds that | 1 | ≤ 1, and there is  < 1 such that | 2 | ≤ , so that As a consequence, the unconditional  2 stability of the scheme is proved.
Remark 1.There are two especial cases of the DPL model that could be singled out.When   = 0, the DPL model reduces to the single-phase-lag (SPL) model [2,22], and (2) reduces to the Cattaneo-Vernotte (CV) model [23,24].In the case where the two phase-lags are exactly the same,   =   , there would be an instantaneous response between the temperature gradient and the heat flux, and (1) would be equivalent to the classical Fourier law, except for a shift in time.However, we note that, regarding the numerical scheme (4) and its stability properties, both especial cases are included in the general stability analysis presented above.

Particular
Case   = 0.In this special case the method is a two-level scheme: so that the amplification matrix reduces to the amplification factor: where with Thus, introducing the parameter  = 1 − cos  ≥ 0, it is not difficult to get to the expressions and from this Therefore, the unconditional  2 stability of the scheme in this special case is also proved.

Conclusions
The use of non-Fourier models of heat conduction in engineering problems requires the use of efficient methods to compute numerical solutions, and a basic condition for these numerical methods to be reliably employed is to be confident that they present good stability properties.The numerical solution procedure for the DPL equation proposed by McDonough et al. [18] provided the base for an efficient numerical solver for higher dimension problems, but its stability properties were not firmly established.
In this note, by analyzing the amplification matrix of the scheme, we have provided a detailed proof showing that the method is indeed unconditionally stable, including its possible application to the particular case when   = 0.
As stressed in [3], the use of appropriate boundary conditions is a key issue in the modelling of non-Fourier heat conduction problems.The type of stability analysis presented in this work is valid for different type of boundary conditions, when they are consistently incorporated into the numerical scheme, and it can also provide correct results in a wide range of situations for nonperiodic boundary conditions [25,Ch. 10].In micro-and nanoscale problems, mixed-type Robin boundary conditions are needed to account for temperature jump on boundaries [26,27], and the incorporation of these types of boundary conditions might require a particular or complementary stability analysis.