Dual-phase lag equation. Stability conditions of a numerical algorithm based on the explicit scheme of the finite difference method
Journal of Applied Mathematics and Computational Mechanics
DUAL-PHASE LAG EQUATION. STABILITY CONDITIONS OF A NUMERICAL ALGORITHM BASED ON THE EXPLICIT SCHEME OF THE FINITE DIFFERENCE METHOD
Ewa Majchrzak 1, Bohdan Mochnacki 2
1 Silesian University of Technology
2 University of Occupational Safety Management in Katowice
Abstract. The dual-phase lag equation (DPLE) is considered. This equation belongs to the group of hyperbolic PDE, contains a second order time derivative and higher order mixed derivative in both time and space. From the engineer’s point of view, the DPLE results from the generalized form of the Fourier law. It is applied as a mathematical model of thermal processes proceeding in the micro-scale and also in the case of bio-heat transfer problem analysis. At the stage of numerical computations the different approximate methods of the PDE solving can be used. In this paper, the authors present the considerations concerning the stability conditions of the explicit scheme of finite difference method (FDM). The appropriate conditions have been found using the von Neumann analysis. In the final part of the paper, the results of testing computations are shown.
Keywords: heat conduction, dual-phase lag equation, finite difference method, stability conditions of FDM explicit scheme
Heat conduction processes proceeding in the domain of a solid body can be described by the Fourier-type equation. As is well known, the mathematical form of this equation results from the assumption of the infinite velocity of thermal wave propagation. In the case of micro-scale heat transfer, and also in the case of materials with a specific internal structure (e.g. biological tissue), the Fourier equation should be modified. To take into account the delay effect of the local and temporary heat flux with respect to the temperature gradient, the so-called relaxation time τq is introduced, and then one obtains the Cattaneo-Vernotte equation . If additionally one assumes the non-zero value of the temperature gradient lag time (the thermalization time τT) then the heat conduction is described by the dual-phase lag equation (e.g. [2, 3]) which is a subject of discussion presented in this paper. The considerations concerning the macro-scale heat conduction do not require the inclusion of lag times (for example, in the case of gold τq = 8.5 ps (1 ps = 10–12 s), τT = 90 ps ), while the duration of the typical heating/cooling process is measured in seconds, at least. However, considering some of the unusual materials (e.g. soft tissue) the lag times should be taken into account despite the fact that one considers the problem in the macro-scale. For example, the experimental investigation made by Antaki  showed that the value of relaxation time is about 2 s for processed meat. Similar problems are discussed in .
The primary aim of this publication is to establish the stability conditions for DPLE (in the case when at the stage of numerical modeling the explicit scheme of the FDM is used). In the opinion of the authors, the conditions reported in the literature are either incomplete or contain errors. Often a time step is chosen intuitively (through trial and error). Here the appropriate conditions have been found using the von Neumann analysis , the testing computations are also performed.
2. The FDM equation and stability conditions
To simplify the mathematical manipulations the 1D task is analyzed. So, the dual phase lag equation (the hyperbolic PDE) in the form
is considered. Here a is the thermal diffusivity, τq and τT are the relaxation and thermalization times. The equation (1) is supplemented by the appropriate boundary and initial conditions.
The numerical solution of the problem discussed can be obtained using the explicit scheme of the FDM. Let us consider the 1D differential mesh being the Cartesian product of the geometric and temporal meshes. Both geometric h and time ∆t mesh steps are assumed to be the constant values.
Below, the FDM equation for the internal nodes (i = 1, 2, ..., n−1) will be presented. The successive differential operators are approximated in the following way:
The approximation error carried by at every node of space i and time f is assumed to have a wave form with the wave number denoted by k and the amplitude by δ :
As time progresses, to assure convergence, the amplitude of the approximation error must be less than unity, i.e. .
Let us introduce the formula (7) into the FDM equation (6)
or (after the simplifications)
After using the Euler formulas one obtains
Let us denote
According to , the absolute values of the roots of equation (12) will be less than 1 when
The solution of the first inequality is of the form
For the worst case (from the viewpoint of stability) one should assume that and then
In the case of the second inequality, one obtains
Obtaing the solutions of the inequalities (14) and (15) is very tedious, of course, therefore only the final results are presented above. Additionally, one can see that the second condition is more ‘demanding’ and in this situation the inequality (18) determines the value of the critical time step.
3. Results of computations
The thin metal film (gold) of thickness 100 nm is considered, The initial temperature Tp = 20°C, the initial heating rate vp = 0. At the moment t = 0 the boundary temperature for x = 0 increased to TB = 50°C, while the temperature for x = G is still equal to Tp (in other words the Dirichlet conditions are assumed). The course of material heating is analyzed. Thermophysical parameters of the film are the following: λ = 315 W/(mK), c = 2.4897 MJ/(m3 K) (this means a = 0.0001265 m/s), τq = 8.5 ps, τT = 90 ps. The mesh step h = 5 nm. The critical time step resulting from inequality (18) is close to ∆tcr = 0.0093 ps. The presented solution in Figure 1 corresponds to ∆t = 0.005 ps and shows the heating curves at the points x = 10, 50, 80 nm. In Figure 2, the unstable solution (temperature profile in the film domain for time t = 1.2 ps) can be observed. The time interval ∆t = 0.01ps is bigger than the critical one and the effect of instability is clearly visible.
Fig. 1. Example of stable solution
Fig. 2. Example of unstable solution
Fig. 3. Stable solution (no-flux condition)
Fig. 4. Unstable solution (no-flux condition)
The solutions of the similar problem are presented in Figures 3 and 4. In the place of the Dirichlet condition for x = G, the no-flux condition (∂T/∂x = 0) has been taken into account.
In the paper the stability conditions for the explicit FDM scheme in the case of the DPL equation are discussed. In particular, the considerations resulting from the von Neumann analysis application are presented. Such an approach leads to two conditions (inequalities), but the fulfillment of the second condition (18) also provides the fulfillment of the first condition (17). It should be also pointed out that the numerical values of the left hand sides of formulas (17) and (18) are very close, but the value of (17) is slightly larger. The condition (18) and its modified form for the case of the control volume method are presented in [9-11]. They have been obtained by transferring a known stability condition for the parabolic equations on the hyperbolic ones (numerical factor C1 - see eq. (5) must be positive). Such an approach has been rather intuitive. The condition (17) was not reported.
The interesting case is shown in Figure 5. The heating curves at the points x = 10, 50, 80 nm have been found for ∆t = 0.009336 ps and this time interval gives a positive value of inequality (17) but a negative one of inequality (18). Despite this fact, the solution obtained is physically correct. This is probably caused from the fact that the critical time step was only minimally exceeded.
Further research will concern the stability conditions for the explicit scheme of the generalized finite difference method  (the DPLE will be considered, of course).
Fig. 5. Solution (temperature profiles) close to the critical time step. The Dirichlet boundary conditions
This work is supported by the project No. 2015/19/B/ST8/01101 and sponsored by the National Science Centre (Poland).
 Cattaneo M.C., A form of heat conduction equation which eliminates the paradox of instantaneous propagation, C.R. Acad. Sci. I - Math. 1958, 247, 431-433.
 Tzou D.Y., Macro- to Microscale Heat Transfer: The Lagging Behavior, John Wiley & Sons, Ltd. 2015.
 Zhang Z.M., Nano/Microscale Heat Transfer, McGraw-Hill, New York 2007.
 Dai W., Nassar R., A domain decomposition method for solving three-dimensional heat transport equations in a double-layered thin film with microscale thickness, Numerical Heat Transfer 2000, Part A, 38, 243-255.
 Antaki P.J., New interpretation of non-Fourier heat conduction in processed meat, ASME J. Heat Transfer 2005,127, 189-193.
 Roetzel W., Putra N., Das S.K., Experiment and analysis for non-Fourier conduction in materials with non-homogeneous inner structure, Int. J. Therm. Sci. 2003, 42, 541-552.
 McDonough J.M., Kunadian I., Kumar R.R., Yang T., An alternative discretization and solution procedure for the dual phase-lag equation, Journal of Computational Physics 2006, 219, 163-171.
 Young D.M., Iterative Solution of Large Linear Systems, Academic Press, New York 1971.
 Majchrzak E., Mochnacki B., Sensitivity analysis of transient temperature field in microdomains with respect to the dual-phase-lag model parameters, Journal for Multiscale Computational Engineering 2014, 12(1), 65-77.
 Mochnacki B., Ciesielski M., Micro-scale heat transfer. Algorithm basing on the Control Volume Method and the identification af relaxation and thermalization times using the search method, Computer Methods in Material Science 2015, 15, 2, 353-361.
 Ciesielski M., Mochnacki B., Application of the control volume method using the Voronoi polygons for numerical modeling of bio-heat transfer problems, Journal of Theoretical and Applied Mechanics 2014, 52, 4, 927-935.
 Mochnacki B., Majchrzak E., Numerical modeling of casting solidification using generalized finite difference method, THERMEC 2009, PTS 1-4 Book Series: Materials Science Forum 2010, 638-642, 2676-2681.