Analysis and Modelling of Non-Fourier Heat Behavior Using the Wavelet Finite Element Method

Non-Fourier heat behavior is an important issue for film material. The phenomenon is usually observed in some laser induced thermal responses. In this paper, the non-Fourier heat conduction problems with temperature and thermal flux relaxations are investigated based on the wavelet finite element method and solved by the central difference scheme for one- and two-dimensional media. The Cattaneo–Vernotte model and the Dual-Phase-Lagging model are used for finite element formulation, and a new wavelet finite element solving formulation is proposed to address the memory requirement problem. Compared with the current methodologies for the Cattaneo–Vernotte model and the Dual-Phase-Lagging model, the present model is a direct one which describe the thermal behavior by one equation about temperature. Compared with the wavelet method proposed by Xiang et al., the developed method can be used for arbitrary shapes. In order to address the efficient computation problems for the Dual-Phase-Lagging model, a novel iteration updating methodology is also proposed. The proposed iteration algorithms on time avoids the use the global stiffness matrix, which allows the efficient calculation for title issue. Numerical calculations have been conducted in the manner of comparisons with the classical finite element method and spectral finite element method. The comparisons from accuracy, efficiency, flexibility, and applicability validate the developed method to be an effective and alternative tool for material thermal analysis.


Introduction
The Fourier law, one of the most important laws to explain the behaviors in heat transport, seems to be ineffective in some researches when high-intense and ultra-short lasers are used as the excitation for microscale heat transports. In these experimental observations of some film materials, the sharp wave fronts responsible for temperature overshooting are hard to be interpreted by the classical model. To address this problem, some new models are proposed, like the Cattaneo-Vernotte (CV) model and the Dual-Phase-Lagging (DPL) model. The macroscopic thermal wave model was firstly postulated by Cattaneo and Vernotte in 1958. The modification leads to a hyperbolic heat condition equation and suggests describing the heat transport by wave with finite speed. However, the two-step model and the pure phonon field model proposed later suggest that the microscale thermal behavior neither follow the pattern given by the CV model nor Fourier diffusion model. To fill the gap between microscopic to macroscopic theories, the DPL model was proposed by Joseph [1] and Tzou [2,3] according to two time constants in the thermal evolution equation. The DPL model aims to remove the precedence assumption made in the CV model. It allows either the temperature gradient to precede the heat

Cattaneo-Vernotte Model (CV Model)
In classical heat conduction investigations, the diffusion of heat is characterized by the empirical law (Fourier law of heat conduction), which postulates that the heat flux is directly proportional to the temperature gradient as: q(r, t) = −k∇T(r, t) where q is the heat flux, r is the position vector, t is the physical time, k is the conductivity for thermal medium, and T is the temperature. Since the diffusion equation is parabolic by nature, it is easy to know from the wave motion view that Equation (1) implies an infinite speed of the propagation of thermal wave, which indicates that the local change in heat flux q can lead to an instantaneous perturbation in temperature field T. It has been verified that the conclusion is incompatible with experiments. Based on the two-fluid model, Tisza [38] predicted the existence of thermal wave, which was detected by Peshkov [39] as the "second sound". Associated with the development of material processing via pulsed sources and the requirement of laser induced guide wave in structural health monitoring, the classical Fourier's law was shown to be inadequate in modelling the high frequency response. The problems mentioned above triggered many attempts to improve the classical model, the most famous one appears to be the CV model, in which the thermal "inertia" is taken into account [7,8]: Equation (2) is also the first order approximation to the single-phase-lag constitutive relation: relaxation time is defined by τ 0 . The physical meaning of τ 0 can be interpreted as the natural result of communication time in molecules collisions in material, which is further formulated as: where α = k/ρc is the thermal diffusivity, ρ and c are the density and the specific heat of the material. The introduction of relaxation time allows the time lag between heat flux and the change of temperature, and hence the thermal wave propagation can be described by this model. Equation (3) and the energy equation: ρc ∂ ∂t T(r, t) + ∇ · q(r, t) = Q(r, t) where Q(r, t) depicts the heat source, yield the delay heat equation: However, the initial value problems for Equation (6) are ill-posed. Therefore, the single-phase-lag constitutive relation (Equation (3)) cannot be considered as sensible physical one [7,8]. In the present references, the Jeffreys-type constitutive relation [1] is usually used to approximate the thermal behavior: τ 0 ∂ 2 ∂t 2 T(r, t) + ∂ ∂t T(r, t) − α∆T(r, t) = 1 ρc Q(r, t) + τ 0 ρc ∂ ∂t Q(r, t) The initial value problems of Equation (7) are well-posed. It should be emphasized that the Jeffreys-type [9,10] constitutive relation Equation (7) cannot be considered as a real approximation of the single-phase-lag relation leading to ill-posed problems.

Dual-Phase-Lagging Model (DPL Model)
It has been confirmed by experimental data that the CV model performs better than the classical Fourier law in numerical prediction. The CV model, however, may obtain some predictions which cannot be supported by experiments. A comprehensive study shows that the CV model has only taken account of the fast-transient effects, but not the micro structural interactions. These two effects can be reasonably represented by the DPL between q and ∇T: where τ T is the delay time caused by the micro-structural interactions such as phonon-electron interaction or phonon scattering and is called the phase-lag of the temperature gradient. Similarly, we can obtain the corresponding delay heat equation by Equations (5) and (8): Likewise, the initial value problems of Equation (9) are ill-posed. The Jeffreys-type constitutive relation of the DPL model is given by: The initial value problems of Equation (10) are also well-posed, thus the Jeffreys-type constitutive relation cannot be considered as a strict description of the DPL model. The higher-order approximations of CV and DPL model were also considered in literatures see in Prof. Rukolaine's and Prof. Chirita's works [7,[40][41][42][43].

The Dimensionless Formulation
Consider that the parameters involved in calculations are very extreme, which generate the difficulty in simulations. The CV and DPL models are usually transformed into the corresponding normalized forms. Firstly, the excitation is nondimensionalized following the formulation introduced in [23]. Traditionally, the Gaussian profile is used to simulate the light intensity of laser pulses: where R is the reflectivity of irradiated surface, I 0 is the output intensity of laser, t p is the full-width-at-half-maximum of pulse. For convenience in the subsequent analysis, the following dimensionless parameters are introduced following the definition in [23]: Time parameters : Temperature parameter : Heat flux parameter : Heat source parameter : where T 0 is the reference temperature. Based on the above dimensionless parameters, Equation (7) for the CV model is rewritten as the dimensionless form: and the dimensionless DPL model of Equation (10) is rewritten as: It is easy to know the CV model can be obtained by the DPL model by setting λ T = 0, namely τ T = 0, in Equation (18). Thus, only the derivation of Equation (18) is presented here. Firstly, substituting the length parameters in Equation (10) by dimensionless forms (either for variables or operators), we get (here R is used to define the dimensionless r): Thereafter, substituting the time parameters in Equation (19) by dimensionless forms: then replace the temperature parameter by dimensionless parameters, considering the relation that α = k/ρc, we get: Lastly, replace the heat source parameter by dimensionless parameters: which is simply denoted in the form given in Equation (18). In above equations, the dimensionless heat source is given by (combine Equations (11) and (16)):

Wavelet Interpolating/Shape Function
In above sections, the basic CV and DPL models are presented in the form of partial differential equations (PDE). In the past decades, many excellent approaches have been developed to obtain the close-form solution of these PDEs, the typical works are referred to Tzou's work and Wang's work [2,11]. The numerical methods, especially the finite element method, however, is more flexural for complex boundary conditions and modelling. The PDEs are then transformed into the formulations which can be used in the FEM.
Similar as the classical finite element method, the region Ω is firstly meshed in terms of a set of nonoverlapping sub-domain Ω e , and each sub-domain is mapped to a unit interval considering the dimension of the problem analyzed. In the unit interval, some mth-order j scale B-spline wavelets on interval φ j m,k (ξ) (BSWI), are used to construct wavelet finite element formulations for title problem. According to the mth-order 0 scale B-spline functions and the corresponding wavelets given by Goswami [44], the j scale mth-order BSWI, which is simply denoted as BSWImj, can be defined. The support of the inner B-spline occupies m segments: At any scale j, the discretization step is 1/2 j . Thus, in order to have at least one inner B-spline function, the following condition should be satisfied: Let j 0 be the initial scale, then for each j ≥ j 0 , The scaling functions φ j m,k (ξ) (m ≥ 2) can be derived by the following formulas [44]: where the function (x) + max(0, x), and x j k , x j k+1 , . . . , x j k+m x is the mth-order divided difference of (x − ξ) m−1 + with respect to x. Referring to Equation (26), we can derive any φ j m,k (ξ) from φ 0 m,k (ξ). Based on Equation (24), the BSWI4 0 (m = 4, j = 0) functions are calculated [44], and then the expressions of BSWI4 3 (m = 4, j = 3) scaling functions are obtained and used as the main interpolating function and shape function in WFEM for example. Restricted by space, we cannot present the BSWI4 3 scaling functions in detail. For two-dimensional case, we further define the horizontal and vertical interpolating vectors based on Equation (27): where ξ, η belong to the interval [0, 1], which depict the normalized x and y coordinates, respectively. The two-dimensional interpolating function is formulated based on the Kronecker product (⊗) between the two vectors in Equation (28), namely Φ = φ ξ ⊗ φ η . Figure 1 presents some typical BSWI4 3 functions in two-dimensional analysis.  In the frame of finite element method (FEM), the unknown continuous temperature field function T(ξ, η, t) can be interpolated in elemental region as: (29) where N is the interpolating function, and T e is the nodal temperature in an element. Since there are more than one node in an element, interpolating function and nodal temperature are both written in matrix form. In this work, the BSWI4 3 function is selected as the interpolating function N. Since the physical filed is recorded in terms of wavelet coefficients in wavelet interpolations, an additional transform matrixT is required to transform the wavelet coefficients into the physical domain. Thereafter, the interpolating function N yields:

WFEM Formulation
The PDE given for the CV model is transformed into the WFEM formula in this section by aid of the trial function. Equation (17) is the strong form of CV model, where the variables (X, Y, λ) for temperature and heat source are abbreviated for the sake of convenience. The requirement of the two-order continuity of Θ, namely the component ∆Θ, compounds the difficulty of trial function selection in WFEM. Therefore, the weak form is usually used. By dotting Equation (17) with the trial function ϑ, and integrating it by parts over the region of interest Ω, one can get: Based on the Hamilton's principle, the weak form of CV model can be obtained in a matrix form, which yields: where the matrixes are defined by: The symbol e defines the total number of finite elements used in modelling, i and j are the indexes of element on different directions, w i and w j are the corresponding weights of Gaussian integrations, and matrix J is the Jacobi matrix. The calculation methodology for these parameters are referred to the basic theory of FEM [45], thus they are not presented in details. Since the structure of Equation (32) is same with the typical wave propagation equation or dynamic responses in elastic problems, it can be predicted that the temperature change propagates in the wave-like mode. The same process of the CV weak form is used to generate the corresponding weak form for the DPL model, where M, C, K, and G have been defined in Equations (33)- (36): The difference between the DPL and the CV model is mainly focused on H, a matrix who plays the role of damping but derived from K. The ratio of relaxation time τ T /τ 0 determines the properties of the DPL model: Equations (32) and (38) present the basic solving formulation of WFEM for the CV and DPL PDEs, however, these basic formulations are only suitable for the small-scale computation, namely the degrees of freedoms (Dof ) are strictly restricted because of limitation of computer memory. To address this problem, a special solving methodology is proposed in this work.

Solving Methodology
The direct solving methodologies for Equations (32) and (38) are well-established in numerical analysis, such as the mode superposition scheme and the central difference time integration scheme. These methodologies, however, are restricted in small Dof issue due to the limitation of the hard memory of computer. In order to interpret this problem, we can consider 1000 Dofs structure to be inspected by CV or DPL model. The Dofs considered here is still a relatively small scale for calculation. By aid of the weak forms, the 1000 × 1000 matrixes M, C, H, and K, and the 1000 × 1 vector G can be obtained. Firstly, one can try to solve Equations (32) and (38) by mode superposition method. In this scheme, the inverse of a 1000 × 1000 matrix should be calculated in advance to get the 1st-1000th thermal modal shapes. Thereafter, the mode superposition can be conducted. Totally, the inverse of a 1000 × 1000 matrix cannot be directly calculated in a common computer due to the problem of efficiency, the advanced numerical method like the Lancozs method should be utilized. It has been proven that the mode superposition method is inefficient for large Dofs problem. For wave propagation problem, especially the thermal wave problem, the direct iteration on time like central difference time integration scheme is more efficient. Unluckily, a very large amount of Dofs should be considered to model and simulate the thermal wave problem due to the wide band of excitation and the tiny relaxion time. The total number of Dof usually achieves 10,000 for a reasonable response and avoiding the numerical oscillation, and the matrixes are all larger than 10,000 × 10,000. The central difference time integration scheme avoids the calculation of inverse matrix, however, "out of memory" still threatens the calculation. In addition, the large Dofs used in modelling to capture the high frequency wave characteristics further compounds "out of memory" problem due to the least stable time interval, which is determined by the Von Neumann conditioning number of the inspected model. For above reasons, we selected the central difference time integration scheme, however, with a special solving methodology to address the "out of memory problem" for title issue. It should be emphasized 10,000 even larger is not problem for some commercial software like ANSYS. We take this example here to say that direct calculation of inverse of a very large matrix is impossible in practice. To calculate the inverse of a large matrix is time-consuming, which in further compounded by the iteration on time.
(1) For the CV model. To solve Equation (32)  on the global matrix level. Assume the equations as described in Equations (40) and (41) according to the central difference time integration scheme: Substituting Equations (40) and (41) for the corresponding derivatives in Equation (32): where ∆λ is the step between the neighbor integration slice in time domain. To address the "out of memory problem", Equation (42) is further rewritten as Equation (43) for the CV model: Note matrixes M and C are in a diagonal form, the inverse of matrixes M 0 , M 1 , and M 2 can be easily obtained following the example: Matrixes M 0 , M 1 , and M 2 are stored in vector forms due to their diagonal property. Thereafter, we now focus on the componentF, which will be decomposed into the elemental level for calculation to avoid global multiplication. The specific algorithm is described in Table 1 referring to [46].
In the above algorithm, calculations are mainly conducted on the elemental level, and thus the "out of memory problem" is addressed.
(2) For the DPL model. The difficulty of solving Equation (38) based on the algorithm presented in Table 1 is generated by the non-diagonal property of matrix H. To address this issue, we proposed the following predict-update format. According to the central difference time integration scheme, Equation (38) is rewritten as: (45) where the matrix H is replaced by λ T K. Due to the non-diagonal property of K, the inverse of matrix M 0 cannot be obtained following the example presented in Equation (43). To address this problem, we rewrite the non-diagonal part to the right of equation as: (46) Thereafter, the inverse of M 0 can be efficiently computed. However, the appearance of the variable to be calculated, namely Θ λ+∆λ , makes it impossible for solving. In the predict-update iteration method, the value of the Θ λ+∆λ on the right-hand side is assigned to be Θ λ in the predicting phase: Equation (46) is modified as (the predicting phase): 1 The solving methodology of Equation (49) is similar with the algorithm shown in Table 1, and then one can get the predicted value Θ p λ+∆λ . Substitute Θ p λ+∆λ for the Θ λ+∆λ on the right-hand side of Equation (49), the updating phase is obtained: (51) Solving Equation (51), the updated T u t+∆t is obtained, and then substituted into Equation (46) for iteration until the convergence.

Definition of Boudary Condition and Initial Condition
Boundary condition and initial condition are important for thermal analysis. The thermal boundary conditions can be clustered into three types: (1) Dirichlet (the first type) boundary condition: the temperatures of some areas are given. (2) Neumann (the second type) boundary condition: the heat fluxes of some areas are given. (3) The mixture of (1) and (2) types. To implement the boundary condition on PDEs may be problematic for these quantities in parallel. However, one can define them in the above algorithms easily. Taking the first type and the second type boundary condition for example: (1) the first type boundary condition, in which temperature is fixed in special area. Let us pay attention to Table 1, in the last command in the loop, we calculated the temperature for time step 1, which is further used as the initial condition for time step 2. Therefore, one can restrain (i.e. set to be the given value) the corresponding temperature for restrained area to satisfy the first type boundary condition. The similar process is repeated in every time step and thus the boundary condition can be fulfilled. (2) the second type boundary condition, in which heat flux is fixed in special area. It is worth to point out that the accurate definition of the second type boundary condition in the present algorithm is impossible. It seems that we can implemented the first order derivative of temperature by Equation (40), by which the heat flux can be restricted. However, the heat flux in the non-Fourier problems, is no longer directly proportional to the temperature gradient, especially in heterogeneous materials or for low-temperature phenomena [47][48][49][50][51]. Of course, one may consider a higher order modification based on Equations (40) and (41), the first and the second order derivatives of temperature. But the method is still not accurate enough. For this reason, heat flux is usually treated as an independent state variable in calculations. Recently, an implicit scheme are presented by Reith and Kovács et al. [52], which successfully describes numerous beyond Fourier experimental findings. More details of the boundary problem can be referred to this work. The restrain of initial condition is to set the initial value to Θ −1 and Θ 0 , and thus the initial value of temperature and heat flux can be controlled.

Stability Conditions of Central Difference Time Integration
A disadvantage of the central differences method can be its conditional stability, which requires that the length of the time step ∆λ be smaller than some critical value that is closely related to the dynamic properties of the discretized system: where ω n is the shortest period of eigenvalue of the discrete system. The stable value of ∆λ the finite element system is determined by the maximum of ω n in system. In practice, the value of ∆λ is usually selected by trial-and-error method, the initial value of ∆λ is set as: ∆λ ≤ ∆s min c max (53) where ∆s min is the minimum characteristic length of element. For one-dimensional case, it is the length of element, and it is the 2-norm (alternative) of all elemental edges for two-dimensional case. When the trial value induces an unstable solution, reduce the ∆λ and repeat the calculation.

Numerical Results and Discussions
Utilizing the CV/DPL models and the corresponding solving formats, numerical computations are performed to display the lagging thermal behavior in various media under pulse-laser heating in the form defined in above sections. Mentioned that the situation λ T = 0 degenerates the DPL into the CV model, the CV model is deemed as a special case of the DPL model and thus different models are not emphasized in the analysis. In the following parts, discussions are organized by the following logic: The boundary conditions of the following cases are all defined as the Neumann boundary heat flux is 0. For one-dimensional case, the boundary is set on the two tips of region. In two-dimensional case, the Neumann boundary is restrained on the boundary of region. For all cases, the 0 initial condition is applied, namely temperature change and the heat flux are all 0 at beginning. The physical parameters before dimensionless process are given for possible comparison: (1) the material with the thermal parameters α = 0.2301 × 10 -4 m 2 /s, τ 0 = 0.1720 ps, (2) geometrical parameters, x = 5 nm (equivalent to 1.257 in dimensionless domain) for one-dimensional case, radius r = 6 nm (equivalent to 1.508 in dimensionless domain) for two-dimensional case. 3) excitation parameter, τ p = 100 fs (equivalent to λ p = 0.2907), the reflectivity of irradiated surface is simply assumed as R = 0. 4) Temporal parameter, for one-dimensional case, t = 3 ps (equivalent to λ = 8.7208), divided into 10,000 time steps; for two-dimensional case t = 1.5 ps (equivalent to λ = 4.3604), divided into 10,000 time steps.

Convergence and Accuracy
The performances of the presented method on convergence and accuracy are validated by the comparison with SFEM and FEM on one-dimensional problem. As the dimensionless parameters are used, they are not emphasized here. The overview of the WFEM performance on one-dimensional thermal wave behavior can be found in Figure 2, where the thermal behavior on the total field are presented for different λ T . However, it is observed that the differences contained in the sub-figures among different λ T is not evident. Thus, Figures 3-5 are further presented for comparison. In this numerical example, the heat source is implemented on X = 0. The total length of the inspected region is 1.257. For different cases, the WFEM shows a good applicability. Covered by the strong thermal impulse, the differences induced by the variation of λ T are not clear on the interval X ⊂ [0, 0.5]. However, on the interval X > 0.5, the dimensionless temperature change becomes smoother with the increase of λ T , which means the diffusive behavior becomes more dominant and the wave behavior becomes weaker.

Convergence and Accuracy
The performances of the presented method on convergence and accuracy are validated by the comparison with SFEM and FEM on one-dimensional problem. As the dimensionless parameters are used, they are not emphasized here. The overview of the WFEM performance on one-dimensional thermal wave behavior can be found in Figure 2 A detailed comparison on accuracy and efficiency is given in Figure 3   A detailed comparison on accuracy and efficiency is given in Figure 3 based on the same numerical scenario λ T = 0 used above. The compared methodologies are FEM, the classical method, SFEM, one of the best methods for dynamic analysis, and the presented WFEM. The FEM is conducted using the classical elements with Dof = 11, 51, and 101, and the compared SFEM and WFEM methods are conducted with the Dof = 11. The temperature changes on excited point (X = 0), middle of the media (X = 0.6283) and the tip of one-dimensional region (X = 1.257) are calculated by above methods. In Figure 3, the poor convergency of FEM is totally performed on this issue. The calculations based on Dof = 11 of FEM show the obvious numerical oscillation on different positions. At point X = 0 (left picture of Figure 3), the numerical oscillation of FEM results has essentially changed the real waveform. The dense grid can suppress this phenomenon, the FEM result of Dof = 101 is still beset by the oscillation (seen in the left picture in Figure 4). At point X = 0.6283 (middle picture of Figure 3), the wave behaviors are not clearly performed by FEM. Either the behavior the direct wave (near λ = 2), or the arrival time of the primary reflection from the tip (near λ = 4) cannot be accurately described by the FEM. At point X = 1.257 (right picture of Figure 3), even a false waveform can be observed near λ = 1.5 in the responses obtained by FEM. Thus, the classical FEM cannot qualify the analysis of the title problem. For this reason, the FEM results for two-dimensional cases are not further presented. Likewise, the WFEM (Dof = 11) and SFEM (Dof = 11) are introduced into analysis. The comparisons on accuracy and efficiency among FEM, SFEM and WFEM are given in Figure 4, where the zoom-out views on some key time intervals are presented for comparison. We can see the SFEM and WFEM achieve the higher accuracy with less Dofs used on the three inspected positions compared with FEM. The agreement between WFEM and SFEM further validates the effectiveness of the developed method. The comparison given in Figure 4 proves the WFEM to be an alternative methodology for thermal wave analysis.  Likewise, the WFEM (Dof = 11) and SFEM (Dof = 11) are introduced into analysis. The comparisons on accuracy and efficiency among FEM, SFEM and WFEM are given in Figure 4, where the zoom-out views on some key time intervals are presented for comparison. We can see the SFEM and WFEM achieve the higher accuracy with less Dof s used on the three inspected positions compared with FEM. The agreement between WFEM and SFEM further validates the effectiveness of the developed method. The comparison given in Figure 4 proves the WFEM to be an alternative methodology for thermal wave analysis. Likewise, the WFEM (Dof = 11) and SFEM (Dof = 11) are introduced into analysis. The comparisons on accuracy and efficiency among FEM, SFEM and WFEM are given in Figure 4, where the zoom-out views on some key time intervals are presented for comparison. We can see the SFEM and WFEM achieve the higher accuracy with less Dofs used on the three inspected positions compared with FEM. The agreement between WFEM and SFEM further validates the effectiveness of the developed method. The comparison given in Figure 4 proves the WFEM to be an alternative methodology for thermal wave analysis.  The conclusion given by Figure 4, however, is not complete as the applicability of WFEM on λ  The conclusion given by Figure 4, however, is not complete as the applicability of WFEM on T λ is not investigated yet. Thereafter, results in Figure 5 are supplied for this reason. Generally, the The conclusion given by Figure 4, however, is not complete as the applicability of WFEM on λ T is not investigated yet. Thereafter, results in Figure 5 are supplied for this reason. Generally, the WFEM performs well for the inspected cases, no numerical oscillation can be observed. At point X = 0 (left picture of Figure 5), the arrivals of the peak move to the right and the primary reflection near λ = 5 is blurred with the increase of λ T . The diffusive behavior becomes more dominant. The phenomenon and tendency can be clearly observed in the responses on X = 0.6283 and 1.257. On the middle point (middle picture of Figure 5), responses for λ T = 0 present to be a typical wave behavior similar with the elastic wave, and the arrival time of the wave is instinct near λ = 1.25. For diffusive (λ T = 0.5) and over-diffusive cases (λ T = 1.5), the arrival time cannot be defined as the thermal behavior is not wave-like. Compared the responses at X = 0.6283 with X = 1.257 (right picture of Figure 5), we can see the temperature disturbance arrives nearly the same time for these two positions, this suggests that in the diffusive and over-diffusive cases, the wave speed is nearly infinite.

Flexibility and Applicability
The convergence and accuracy of the presented method have been verified in the above section. Hereafter, the flexibility and applicability of the presented method are validated by computations and comparisons on two-dimensional grids as shown in Figure 6. The inspected region with the dimensionless radius r = 1.508 is meshed by 432 and 400 WFEM elements as shown in Figure 6a Figure  5), we can see the temperature disturbance arrives nearly the same time for these two positions, this suggests that in the diffusive and over-diffusive cases, the wave speed is nearly infinite.

Flexibility and Applicability
(a) (b) The convergence and accuracy of the presented method have been verified in the above section. Hereafter, the flexibility and applicability of the presented method are validated by computations and comparisons on two-dimensional grids as shown in Figure 6. The inspected region with the dimensionless radius r = 1.508 is meshed by 432 and 400 WFEM elements as shown in Figure 6a,b, respectively. The heat source locates at point (−1.508, 0), namely the left tip of the region. Although the two cases are meshed in totally different fashions, the similar numbers of elements (432 and 400) are used in calculations to keep the comparability of simulated results. Firstly, the results obtained by the grids shown in Figure 6a,b are presented in Figures 7 and 8 for comparison. It is very clear that we can get the same thermal behaviors via both two kinds of meshes. On another aspect, the quality of grid in Figure 6b is worse than that in Figure 6a since the heterogeneity, especially on the left and right ends. The result of which, however, achieves a similar accuracy with the fine mesh. This can verify the flexibility and applicability of the presented method to some extent.
Hereafter, the flexibility and applicability of the presented method are validated by computations and comparisons on two-dimensional grids as shown in Figure 6. The inspected region with the dimensionless radius r = 1.508 is meshed by 432 and 400 WFEM elements as shown in Figure 6a,b, respectively. The heat source locates at point (−1.508, 0), namely the left tip of the region. Although the two cases are meshed in totally different fashions, the similar numbers of elements (432 and 400) are used in calculations to keep the comparability of simulated results. Then, we will pay attention to the inner comparison of Figure 7, which illustrates the different thermal behaviors in the inspected area. Figure 7a presents the wavy behavior of temperature distribution changes in the inspected area. In this case, the pulsed thermal disturbance propagates in the form of wave. In the snapshot λ = 2.1802, we can observe a clear circle wave front, on which the energy of the pulsed thermal disturbance is mainly focused. With the increase of dimensionless time (λ = 4.3604), the energy is slowly absorbed by the media as the result of matrix C in the DPL model.
With the increase of λ T from 0 to 0.1, the sharp wave fronts in Figure 7a are smoothed and the portions of the disturbance are dissipated. Thus, we can see the influenced area of Figure 7b is larger than that of Figure 7a. Due to the increase of λ T , matrix H plays the more important role in the damping term. As the consequence, the temperature changes on area after wave front becomes more evident compared with the wavy behavior. This is called the wavelike behavior because the concepts like wave front, reflection still can be used to describe the behavior. In both the wavy or wavelike behaviors, some hot points are the reflections of thermal wave to be focused, like the tips of wave fronts. When λ T = 0.5 (Figure 7c), all wavy features disappear, the disturbance caused by pulse transports by diffusion completely. The hot points appear on the tips of wave front are smoothed and the temperature peak in this case is always the heating spot (−1.508, 0). If we continue to enlarge λ T to 1.5, the over-diffusive behavior can be observed (as shown in Figure 7d). Compared with the normal diffusion, the larger λ T produces the higher diffusion rate in an early stage. However, a longer time is required to reach the thermal equilibrium.
To illustrate and compare the thermal behavior shown in Figure 7 more clearly, the thermal responses on points (−1.508, 0), (−0.754, 0), and (0, 0) with different λ T are presented in Figure 9. Totally, the changing trends of the responses in two-dimensional media are similar with that in one-dimensional media, however, with longer time to achieve thermal equilibrium. In addition, the temperature change is also smaller compare with that in the one-dimensional case.
fronts. When 0.5 T λ = (Figure 7c), all wavy features disappear, the disturbance caused by pulse transports by diffusion completely. The hot points appear on the tips of wave front are smoothed and the temperature peak in this case is always the heating spot (−1.508, 0). If we continue to enlarge T λ to 1.5, the over-diffusive behavior can be observed (as shown in Figure 7d). Compared with the normal diffusion, the larger T λ produces the higher diffusion rate in an early stage. However, a longer time is required to reach the thermal equilibrium. To illustrate and compare the thermal behavior shown in Figure 7 more clearly, the thermal responses on points (−1.508, 0), (−0.754, 0), and (0, 0) with different T λ are presented in Figure 9.
Totally, the changing trends of the responses in two-dimensional media are similar with that in onedimensional media, however, with longer time to achieve thermal equilibrium. In addition, the temperature change is also smaller compare with that in the one-dimensional case.

Conclusions
In order to describe the thermal behavior of film material, a new WFEM formulation for thermal behaviors in one-and two-dimensional media has been proposed in this work. The hyperbolic heat conduction model has been solved using the central difference scheme on time and wavelet interpolation in space. The proposed algorithm has been tested by comparison with the classical FEM and SFEM on the aspects of accuracy and efficiency. The flexibility and applicability for different mesh grids are validated in a two-dimensional case. The current work provides an alternative tool for the analysis of thermal analysis. It should be mentioned that the presented method cannot apply the heat flux boundary condition accurately due to the FEM formulation used in this work. It is worth to investigate the state variable based WFEM further to address this problem.

Conclusions
In order to describe the thermal behavior of film material, a new WFEM formulation for thermal behaviors in one-and two-dimensional media has been proposed in this work. The hyperbolic heat conduction model has been solved using the central difference scheme on time and wavelet interpolation in space. The proposed algorithm has been tested by comparison with the classical FEM and SFEM on the aspects of accuracy and efficiency. The flexibility and applicability for different mesh grids are validated in a two-dimensional case. The current work provides an alternative tool for the analysis of thermal analysis. It should be mentioned that the presented method cannot apply the heat flux boundary condition accurately due to the FEM formulation used in this work. It is worth to investigate the state variable based WFEM further to address this problem.