Next Article in Journal
Characterization of Ablated Bone and Muscle for Long-Pulsed Laser Ablation in Dry and Wet Conditions
Previous Article in Journal
Non-Precious Electrodes for Practical Alkaline Water Electrolysis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
School of Mechanical Engineering, Xi’an Jiaotong University, Xi’an 710049, China
2
The State Key Laboratory for Manufacturing Systems Engineering, Xi’an Jiaotong University, Xi’an 710049, China
*
Authors to whom correspondence should be addressed.
Materials 2019, 12(8), 1337; https://doi.org/10.3390/ma12081337
Submission received: 1 April 2019 / Revised: 21 April 2019 / Accepted: 23 April 2019 / Published: 24 April 2019

Abstract

:
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.

1. 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 flux vector or the heat flux vector to precede the temperature gradient in the transient process [4]. The model tries to lump the microstructural effect into the delay of response in time. Although it is still a postulated model, some experimental observations given by Tzou and Tang et al. [5,6] have shown the well agreement with it. It should be emphasized that the CV and DPL model used in many works are derived using a Taylor series expansion. In fact, this way of derivation is not compatible with the second law of thermodynamics. Moreover, the DPL equation is a special, linearized version of the Jeffrey equation which is compatible with thermodynamics. In Rukolaine’s work [7,8,9,10], it is found that the parameters appearing in the DPL equation are not independent from each other. Thus, they cannot be arbitrary, otherwise the solutions can be unphysical.
Compared with the parabolic diffusion equation of Fourier model, the CV and DPL models are hyperbolic in nature. As a result, there is a resurgent interest in the solution of the heat conduction equations given by the CV or DPL model, which accounts for the finite propagation velocity of a thermal wave within the inspected media A survey of numerical schemes for the solving of the heat conduction equation can be organized into two categories, namely analytical solutions [4,11] and numerical solutions [12,13,14,15]. Some remarkable analytical methodologies including the Laplace transformation [16], Green’s function [17], and the integral equation method [18] have been widely investigated. The details of these analytical solution methodologies can be referred to Wang’s book [11] for more information. Without doubt, the accuracy and efficiency of analytical solution are unparalleled compared to these properties of numerical methods. However, we cannot overlook the sophisticated mathematical skills and the complex transformation required in analytical methods, which have been the main obstacles for their application in practice. Due to the complexity of hyperbolic equation, especially the complex inspected region, only very few simple cases can be solved analytically. Consequently, more attentions have been drawn by the numerical solutions and methodologies. However, the accurate solution for film material and two-dimensional medias are sometimes not easily obtained [19,20]. The finite difference method appears to be firstly used to analyze this problem by Yeung et al. [21,22]. They introduced a simple and concise finite difference algorithm developed by applying the Godunov scheme on the characteristic equation. After that, Han and Tang et al. extended the work to the two-dimensional media [23]. Dai et al. further developed the finite difference algorithm, based on which they proposed the convergent three-level finite difference scheme [24] and a high order accurate finite difference method for solving the two-dimensional DPL problem [25]. Boltzmann method is also a popular method used to solve the title problem. Wang et al. [26] developed an enhanced Gray model by considering the second-order terms in Taylor expansion in the phonon Boltzmann transport equation. Xu and Wang derived the DPL model using the Boltzmann method, and investigated the oscillation of the microscale heat conduction systemically [27,28]. Beside the above methods, the classical finite element method (FEM) is also usually used. Ai et al. [29,30] constructed a discontinuous FEM model and analyzed the thermal wave propagation in one- and two-dimensional medias using DPL model. Motivated by these works, this article aims to develop the wavelet finite element (WFEM) formulation and the corresponding solving methodologies for non-Fourier heat conduction. The wavelet finite element is a novel element proposed by Xiang et al. [31,32]. In the authors’ former works [33,34,35,36], the effectiveness of WFEM has been verified for dynamic analysis and elastic wave problems. However, the difference between former issues with the thermal problem make it necessary to develop the WFEM. Now the WFEM is seldomly used in thermal analysis, the only available works are developed by Zhao et al. [37]. Zhao’s works show the potential of WFEM for thermal analysis. Considering the non-Fourier heat conduction problems have not been investigated, we focus on the WFEM formulation for CV and DPL models.
This paper is organized as follow: Section 2 presents the basic models used in this work, based on that, Section 3 presents the numerical formulations for solving by WFEM. The numerical results and discussions are given in Section 4 for validation.

2. Problem Descriptions

2.1. 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]:
τ 0 q ( r , t ) t + q ( r , t ) = k T ( r , t )
Equation (2) is also the first order approximation to the single-phase-lag constitutive relation:
q ( r , t + τ 0 ) = k T ( r , t )
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:
τ 0 = α v C V 2 v C V = α τ 0 = k ρ c τ 0
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:
t T ( r , t ) α Δ T ( r , t τ 0 ) = Q ( r , t )
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.

2.2. 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 :
q ( r , t + τ 0 ) = k T ( r , t + τ 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):
t T ( r , t ) α Δ T ( r , t τ 0 + τ T ) = Q ( r , t )
Likewise, the initial value problems of Equation (9) are ill-posed. The Jeffreys-type constitutive relation of the DPL model is given by:
τ 0 2 t 2 T ( r , t ) + t T ( r , t ) α Δ T ( r , t ) α τ T t Δ T ( r , t ) = 1 ρ c Q ( r , t ) + τ 0 ρ c t Q ( r , t )
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].

2.3. 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:
Q = ( 1 R ) I 0 π t p exp ( 1 t 2 t p 2 )
where R is the reflectivity of irradiated surface, I0 is the output intensity of laser, tp 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]:
Length   parameters :   X = x 2 α τ 0 ,   Y = y 2 α τ 0 ,   R c = r c 2 α τ 0 ,
Time   parameters :   λ = t 2 τ 0 ,   λ p = t p 2 τ 0 ,   λ T = τ T 2 τ 0 ,
Temperature   parameter :   Θ = π τ 0 α k ( T T 0 ) ( 1 R ) I 0 ,
Heat   flux   parameter :   φ = τ 0 π ( 1 R ) I 0 q ,
Heat   source   parameter :   ψ = 2 τ 0 π α τ 0 ( 1 R ) I 0 Q ,
where T0 is the reference temperature. Based on the above dimensionless parameters, Equation (7) for the CV model is rewritten as the dimensionless form:
2 λ 2 Θ + 2 λ Θ Δ Θ = 2 ψ + λ ψ
and the dimensionless DPL model of Equation (10) is rewritten as:
2 λ 2 Θ + 2 λ Θ Δ Θ λ T λ Δ Θ = 2 ψ + λ ψ
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):
τ 0 2 t 2 T ( R , t ) + t T ( R , t ) 1 4 τ 0 Δ T ( R , t ) τ T 4 τ 0 t Δ T ( R , t ) = 1 ρ c Q ( R , t ) + τ 0 ρ c t Q ( R , t )
Thereafter, substituting the time parameters in Equation (19) by dimensionless forms:
2 λ 2 T ( R , λ ) + 2 λ T ( R , λ ) Δ T ( R , λ ) λ T λ Δ T ( R , λ ) = 4 τ 0 ρ c Q ( R , λ ) + 4 τ 0 2 ρ c λ Q ( R , λ )
then replace the temperature parameter by dimensionless parameters, considering the relation that α = k/ρc, we get:
2 λ 2 Θ ( R , λ ) + 2 λ Θ ( R , λ ) Δ Θ ( R , λ ) λ T λ Δ Θ ( R , λ ) = 4 τ 0 α π τ 0 ( 1 R ) I 0 Q ( R , λ ) + 2 τ 0 α π τ 0 ( 1 R ) I 0 λ Q ( R , λ )
Lastly, replace the heat source parameter by dimensionless parameters:
2 λ 2 Θ ( R , λ ) + 2 λ Θ ( R , λ ) Δ Θ ( R , λ ) λ T λ Δ Θ ( R , λ ) = 2 ψ ( R , λ ) + λ ψ ( R , λ )
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)):
ψ = α τ 0 λ p exp ( 1 λ 2 λ p 2 )

3. Numerical Model

3.1. 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 ϕ m , k j ( ξ ) (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:
{ 0   boundary :   x m + 1 j = x m + 2 j = = x 0 j = 0 inner   knots :   x k j = k 2 j   k = 0 , 1 , , 2 j 1   boundary :   x 2 j + 1 j = x 2 j + 2 j = = x 2 j + m 1 j = 1
At any scale j, the discretization step is 1/2j. Thus, in order to have at least one inner B-spline function, the following condition should be satisfied:
2 j 2 m 1
Let j0 be the initial scale, then for each j j 0 ,
ϕ m , k j ( ξ ) = { ϕ m , k j 0 ( 2 j j 0 ξ ) ,   k = m + 1 , , 1 ϕ m , 2 j m k j 0 ( 1 2 j j 0 ξ ) ,   k = 2 j m + 1 , , 2 j 1   ϕ m , 0 j 0 ( 2 j j 0 ξ 2 j 0 k ) ,   k = 0 , , 2 j m ( 0   boundary   scaling   functions ) ( 1   boundary   scaling   functions ) ( inner   scaling   functions )
The scaling functions ϕ m , k j ( ξ ) ( m 2 ) can be derived by the following formulas [44]:
ϕ m , k j ( ξ ) = ( x k + m j x k j ) × [ x k j , x k + 1 j , , x k + m j ] x ( x ξ ) + m 1
where the function ( x ) + max ( 0 , x ) , and [ x k j , x k + 1 j , , x k + m j ] x is the mth-order divided difference of ( x ξ ) + m 1 with respect to x. Referring to Equation (26), we can derive any ϕ m , k j ( ξ ) from ϕ m , k 0 ( ξ ) . Based on Equation (24), the BSWI40 (m = 4, j = 0) functions are calculated [44], and then the expressions of BSWI43 (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 BSWI43 scaling functions in detail. For two-dimensional case, we further define the horizontal and vertical interpolating vectors based on Equation (27):
ϕ ξ { ϕ m , m + 1 j ( ξ )   ϕ m , m + 2 j ( ξ ) ϕ m , 2 j 1 j ( ξ ) } , ϕ η { ϕ m , m + 1 j ( η )   ϕ m , m + 2 j ( η ) ϕ m , 2 j 1 j ( η ) }
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 BSWI43 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:
T ( ξ , η , t ) = N T e
where N is the interpolating function, and Te 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 BSWI43 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 matrix T ^ is required to transform the wavelet coefficients into the physical domain. Thereafter, the interpolating function N yields:
Φ T ^ = N
where the transform matrix is T ^ = { ϕ ξ T ( ξ 1 ) ,   ϕ x T ( ξ 2 )     ϕ x T ( ξ n + 1 ) } T { ϕ η T ( η 1 ) ,   ϕ y T ( η 2 )     ϕ y T ( η n + 1 ) } T .

3.2. 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:
Ω ϑ 2 λ 2 Θ d Ω + 2 Ω ϑ λ Θ d Ω + Ω ϑ Θ d Ω = Ω ϑ ( 2 ψ + λ ψ ) d Ω
Based on the Hamilton’s principle, the weak form of CV model can be obtained in a matrix form, which yields:
M 2 λ 2 Θ + C λ Θ + K Θ = G
where the matrixes are defined by:
M = e i n + 1 j n + 1 w i w j N T N det ( J )
C = e i n + 1 j n + 1 2 w i w j N T N det ( J )
K = e i n + 1 j n + 1 w i w j N T N det ( J )
G = e i n + 1 j n + 1 w i w j N T ( 2 ψ + λ ψ ) det ( J )
The symbol e defines the total number of finite elements used in modelling, i and j are the indexes of element on different directions, wi and wj 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):
Ω ϑ 2 λ 2 Θ d Ω + 2 Ω ϑ λ Θ d Ω + Ω ϑ Θ d Ω + λ T Ω ϑ λ Θ d Ω = Ω ϑ ( 2 ψ + λ ψ ) d Ω
M 2 λ 2 Θ + ( C + H ) λ Θ + K Θ = G
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:
H = λ T e i n + 1 j n + 1 w i w j N T N det ( J ) = λ T K
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.

3.3. 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) via the central difference time integration scheme, the main difficulty focuses on the decomposition of large matrix to avoid the direct calculation like K Θ on the global matrix level. Assume the equations as described in Equations (40) and (41) according to the central difference time integration scheme:
Θ λ λ = Θ λ + Δ λ Θ λ Δ λ 2 Δ λ
2 Θ λ λ 2 = Θ λ + Δ λ 2 Θ λ + Θ λ Δ λ Δ λ 2
Substituting Equations (40) and (41) for the corresponding derivatives in Equation (32):
M Θ λ + Δ λ 2 Θ λ + Θ λ Δ λ Δ λ 2 + C Θ λ + Δ λ Θ λ Δ λ 2 Δ λ + K Θ λ = G λ
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:
( 1 Δ λ 2 M + 1 2 Δ λ C ) M 0 Θ λ + Δ λ = G λ K Θ λ F ^ + ( 2 Δ λ 2 M ) M 1 Θ λ ( 1 Δ λ 2 M 1 2 Δ λ C ) M 2 Θ λ Δ λ
Note matrixes M and C are in a diagonal form, the inverse of matrixes M0, M1, and M2 can be easily obtained following the example:
( M 0 ) 1 = d i a g ( 1 M i i 0 )
Matrixes M0, M1, and M2 are stored in vector forms due to their diagonal property. Thereafter, we now focus on the component F ^ , 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:
( 1 Δ λ 2 M + 1 2 Δ λ C + λ T 2 Δ λ K ) M 0 Θ λ + Δ λ = G λ K Θ λ F ^ + ( 2 Δ λ 2 M ) M 1 Θ λ ( 1 Δ λ 2 M 1 2 Δ λ C λ T 2 Δ λ K ) M 2 Θ λ Δ λ
where the matrix H is replaced by λ T K . Due to the non-diagonal property of K, the inverse of matrix M0 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:
( 1 Δ λ 2 M + 1 2 Δ λ C ) M 0 Θ λ + Δ λ = G λ K Θ λ F ^ + λ T 2 Δ λ K Θ λ Δ λ λ T 2 Δ λ K Θ λ + Δ λ + ( 2 Δ λ 2 M ) M 1 Θ λ ( 1 Δ λ 2 M 1 2 Δ λ C ) M 2 Θ λ Δ λ
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 Δ λ 2 M + 1 2 Δ λ C ) M 0 Θ λ + Δ λ = G λ K Θ λ F ^ + ( 2 Δ λ 2 M ) M 1 Θ λ ( 1 Δ λ 2 M 1 2 Δ λ C ) M 2 Θ λ Δ λ
( 1 Δ λ 2 M + 1 2 Δ λ C ) Θ λ + Δ λ p = G λ ( 1 + λ T 2 Δ λ ) K Θ λ + λ T 2 Δ λ K Θ λ Δ λ + ( 2 Δ λ 2 M ) Θ λ ( 1 Δ λ 2 M 1 2 Δ λ C ) Θ λ Δ λ
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:
Θ λ + Δ λ Θ λ + Δ λ p
( 1 Δ λ 2 M + 1 2 Δ λ C ) Θ λ + Δ λ u = G λ K Θ λ + λ T 2 Δ λ K Θ λ Δ λ λ T 2 Δ λ K Θ λ + Δ λ p + ( 2 Δ λ 2 M ) Θ λ ( 1 Δ λ 2 M 1 2 Δ λ C ) Θ λ Δ λ
Solving Equation (51), the updated T t + Δ t u is obtained, and then substituted into Equation (46) for iteration until the convergence.

3.4. 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.

3.5. 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:
Δ λ Δ λ c r i t i c a l = 2 ω n
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
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.

4. 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:
(1)
Validate the convergence and accuracy of the presented WFEM method by comparing with the time domain spectral finite element method (SFEM) proposed by Ostachowicz and Kudela et al. [53], one of the best methods for the dynamic analysis, and the classical FEM. The comparisons about convergence and accuracy are conducted on one-dimensional structures. These methods are all coded by Matlab in the similar program structure. It should be mentioned that although the time consumption can be obtained by “tic, toc” in Matlab and the similar program structure are used, we do not compare the efficiency by time, however, by DOFs used.
(2)
Different behaviors of the inspected systems are performed using the developed model, containing the wavy behavior ( λ T = 0 ), the wavelike behavior ( λ T = 0.1 ), the diffusive behavior ( λ T = 0.5 ), and the over-diffusive behavior ( λ T = 1.5 ). On this aspect, the applicability of the proposed model in different situations can be verified.
(3)
Considering the simplicity of one-dimensional grids, the flexibility and applicability of the presented method are validated by comparisons on two-dimensional grids.
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 m2/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.

4.1. 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, Figure 3, Figure 4 and Figure 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.
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.
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.

4.2. 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,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 Figure 7 and Figure 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.
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.

5. 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.

Author Contributions

Conceptualization, Z.-B.Y. and X.-F.C.; Funding acquisition, Z.-B.Y. and S.-H.T.; Investigation, Z.-B.Y. and S.-H.T.; Methodology, Z.-B.Y. and Z.-K.W.; Supervision, X.-F.C.; Validation, Z.-B.Y. and Z.-K.W.; Writing—original draft, Z.-B.Y.; Writing—review and editing, Z.-B.Y.

Funding

This research is funded by the National Natural Science Foundation of China (Nos. 51875433 and 51605365), the Young Talent fund of University Association for Science and Technology in Shaanxi of China (No. 20170502), and the Key Research and development program of Shanxi province (Nos. 2017ZDCXL-GY-02-01 and 2017ZDCXL-GY-02-02).

Acknowledgments

The first author Z.-B.Y. wants to show his special thanks for his lovely sons S.-H. and S.-K. (their mother seems to be more thankworthy to some extent). It is their generous and punctual sleeping recently that allows Z.-B.Y. to have plenty of time on this work.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Joseph, D.D.; Preziosi, L. Heat waves. Rev. Mod. Phys. 1989, 61, 41. [Google Scholar] [CrossRef]
  2. Tzou, D.Y. Macro-to Microscale Heat Transfer: The Lagging Behavior; John Wiley & Sons: London, UK, 2014. [Google Scholar]
  3. Tzou, D.Y. The generalized lagging response in small-scale and high-rate heating. Int. J. Heat Mass Transf. 1995, 38, 3231–3240. [Google Scholar] [CrossRef]
  4. Xu, M.; Wang, L. Thermal oscillation and resonance in dual-phase-lagging heat conduction. Int. J. Heat Mass Transf. 2002, 45, 1055–1061. [Google Scholar] [CrossRef]
  5. Tang, D.; Araki, N. Wavy, wavelike, diffusive thermal responses of finite rigid slabs to high-speed heating of laser-pulses. Int. J. Heat Mass Transf. 1999, 42, 855–860. [Google Scholar] [CrossRef]
  6. Tzou, D.Y. Experimental support for the lagging behavior in heat propagation. J. Thermophys. Heat Transf. 1995, 9, 686–693. [Google Scholar] [CrossRef]
  7. Rukolaine, S.A. Unphysical effects of the dual-phase-lag model of heat conduction: Higher-order approximations. Int. J. Thermal Sci. 2017, 113, 83–88. [Google Scholar] [CrossRef]
  8. Rukolaine, S.A. Unphysical effects of the dual-phase-lag model of heat conduction. Int. J. Heat Mass Transf. 2014, 78, 58–63. [Google Scholar] [CrossRef]
  9. Rukolaine, S.A.; Samsonov, A.M.J.P.R.E. Local immobilization of particles in mass transfer described by a Jeffreys-type equation. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 2013, 88, 062116. [Google Scholar] [CrossRef]
  10. Rukolaine, S.A.; Samsonov, A.M. A model of diffusion, based on the equation of the Jeffreys type. In Proceedings of the International Conference Days on Diffraction 2013, St. Petersburg, Russia, 27–31 May 2013. [Google Scholar]
  11. Wang, L.; Zhou, X.; Wei, X. Heat Conduction: Mathematical Models and Analytical Solutions; Springer: Berlin, Germany, 2007. [Google Scholar]
  12. Zhang, S.; Dai, W.; Wang, H.; Melnik, R.V.N. A finite difference method for studying thermal deformation in a 3D thin film exposed to ultrashort pulsed lasers. Int. J. Heat Mass Transf. 2008, 51, 1979–1995. [Google Scholar] [CrossRef]
  13. Du, X.; Dai, W.; Wang, P. A finite-difference method for studying thermal deformation in a 3-D microsphere exposed to ultrashort pulsed lasers. Numer. Heat Transf. Part A Appl. 2007, 53, 457–484. [Google Scholar] [CrossRef]
  14. Lu, W.-Q.; Liu, J.; Zeng, Y. Simulation of the thermal wave propagation in biological tissues by the dual reciprocity boundary element method. Eng. Anal. Bound. Elem. 1998, 22, 167–174. [Google Scholar] [CrossRef]
  15. Hosseini-Tehrani, P.; Eslami, M. BEM analysis of thermal and mechanical shock in a two-dimensional finite domain considering coupled thermoelasticity. Eng. Anal. Bound. Elem. 2000, 24, 249–257. [Google Scholar] [CrossRef]
  16. Liu, K.-C.; Cheng, P.-J.; Wang, Y.-N. Analysis of non-Fourier thermal behaviour for multi-layer skin model. Therm. Sci. 2011, 15, 61. [Google Scholar]
  17. Mandelis, A. Diffusion-Wave Fields: Mathematical Methods and Green Functions; Springer: Berlin, Germany, 2013. [Google Scholar]
  18. Wu, C.-Y. Integral equation solution for hyperbolic heat conduction with surface radiation. Int. Commun. Heat Mass Transf. 1988, 15, 365–374. [Google Scholar] [CrossRef]
  19. Tornabene, F.; Fantuzzi, N.; Bacciocchi, M. Linear Static Behavior of Damaged Laminated Composite Plates and Shells. Materials 2017, 10, 811. [Google Scholar] [CrossRef] [PubMed]
  20. Tornabene, F.; Fantuzzi, N.; Bacciocchi, M.; Viola, E.; Reddy, J.N. A Numerical Investigation on the Natural Frequencies of FGM Sandwich Shells with Variable Thickness by the Local Generalized Differential Quadrature Method. Appli. Sci. 2017, 7, 131. [Google Scholar] [CrossRef]
  21. Lam, T.T.; Yeung, W.K. A numerical scheme for non-Fourier heat conduction, part II: two-dimensional problem formulation and verification. Numer. Heat Transf. Part B Fundam. 2002, 41, 543–564. [Google Scholar] [CrossRef]
  22. Yeung, W.K.; Lam, T.T. A numerical scheme for non-Fourier heat conduction, part I: one-dimensional problem formulation and applications. Numer. Heat Transf. 1998, 33, 215–233. [Google Scholar] [CrossRef]
  23. Han, P.; Tang, D.; Zhou, L. Numerical analysis of two-dimensional lagging thermal behavior under short-pulse-laser heating on surface. Int. J. Eng. Sci. 2006, 44, 1510–1519. [Google Scholar] [CrossRef]
  24. Dai, W.; Shen, L.; Nassar, R. A convergent three-level finite difference scheme for solving a dual-phase-lagging heat transport equation in spherical coordinates. Numer. Methods Partial Differ. Equ. An Int. J. 2004, 20, 60–71. [Google Scholar] [CrossRef]
  25. Sun, H.; Du, R.; Dai, W.; Sun, Z. A high order accurate numerical method for solving two-dimensional dual-phase-lagging equation with temperature jump boundary condition in nanoheat conduction. Numer. Methods Partial Differ. Equ. 2015, 31, 1742–1768. [Google Scholar] [CrossRef]
  26. Wang, D.; Qu, Z.; Ma, Y. An enhanced Gray model for nondiffusive heat conduction solved by implicit lattice Boltzmann method. Int. J. Heat Mass Transf. 2016, 94, 411–418. [Google Scholar] [CrossRef]
  27. Xu, M.; Wang, L. Dual-phase-lagging heat conduction based on Boltzmann transport equation. Int. J. Heat Mass Transf. 2005, 48, 5616–5624. [Google Scholar] [CrossRef]
  28. Cheng, L.; Xu, M.; Wang, L. From Boltzmann transport equation to single-phase-lagging heat conduction. Int. J. Heat Mass Transf. 2008, 51, 6018–6023. [Google Scholar] [CrossRef]
  29. Ai, X.; Li, B. A discontinuous finite element method for hyperbolic thermal wave problems. Eng. Comput. 2004, 21, 577–597. [Google Scholar] [CrossRef]
  30. Ai, X.; Li, B. Numerical simulation of thermal wave propagation during laser processing of thin films. J. Electron. Mater. 2005, 34, 583–591. [Google Scholar] [CrossRef]
  31. Xiang, J.; Chen, X.; He, Y.; He, Z. The construction of plane elastomechanics and Mindlin plate elements of B-spline wavelet on the interval. Finite Elem. Anal. Des. 2006, 42, 1269–1280. [Google Scholar] [CrossRef]
  32. Xiang, J.; Chen, X.; Mo, Q.; He, Z. Identification of crack in a rotor system based on wavelet finite element method. Finite Elem. Anal. Des. 2007, 43, 1068–1081. [Google Scholar] [CrossRef]
  33. Yang, Z.; Chen, X.; He, Y.; He, Z.; Zhang, J. The Analysis of Curved Beam Using B-Spline Wavelet on Interval Finite Element Method. Shock Vib. 2014, 2014. [Google Scholar] [CrossRef]
  34. Yang, Z.; Chen, X.; Li, X.; Jiang, Y.; Miao, H.; He, Z. Wave motion analysis in arch structures via wavelet finite element method. J. Sound Vib. 2014, 333, 446–469. [Google Scholar] [CrossRef]
  35. Yang, Z.; Chen, X.; Zhang, X.; He, Z. Free vibration and buckling analysis of plates using B-spline wavelet on the interval Mindlin element. Appl. Math. Model. 2013, 37, 3449–3466. [Google Scholar] [CrossRef]
  36. Zuo, H.; Yang, Z.; Chen, X.; Xie, Y.; Miao, H. Analysis of laminated composite plates using wavelet finite element method and higher-order plate theory. Compos. Struct. 2015, 131, 248–258. [Google Scholar] [CrossRef]
  37. Zhao, B. Temperature-coupled field analysis of LPG tank under fire based on wavelet finite element method. J. Therm. Anal. Calorim. 2014, 117, 413–422. [Google Scholar] [CrossRef]
  38. Tisza, L. The thermal superconductivity of helium II and the statistics of Bose-Einstein. Compt. Rend 1938, 207, 1035–1037. [Google Scholar]
  39. Peshkov, V. The Second Sound in Helium II. J. Phys. 1944, 8, 381–382. [Google Scholar]
  40. Chiriţă, S.; Ciarletta, M.; Tibullo, V. On the thermomechanical consistency of the time differential dual-phase-lag models of heat conduction. Int. J. Heat Mass Transf. 2017, 114, 277–285. [Google Scholar] [CrossRef]
  41. Chiriţă, S. On the time differential dual-phase-lag thermoelastic model. Meccanica 2017, 52, 349–361. [Google Scholar] [CrossRef]
  42. Chiriţă, S.; Ciarletta, M.; Tibullo, V. Qualitative properties of solutions in the time differential dual-phase-lag model of heat conduction. Appl. Math. Model. 2017, 50, 380–393. [Google Scholar] [CrossRef]
  43. Chiriţă, S.; Ciarletta, M.; Tibullo, V. On the wave propagation in the time differential dual-phase-lag thermoelastic model. Proc. R. Soc. A Math. Phys. Eng. Sci. 2015, 471, 20150400. [Google Scholar] [CrossRef]
  44. Goswami, J.C.; Chan, A.K.; Chui, C.K. On Solving First-Kind Integral-Equations Using Wavelets on a Bounded Interval. IEEE Trans. Antennas Propag. 1995, 43, 614–622. [Google Scholar] [CrossRef]
  45. Zienkiewicz, O.C.; Taylor, R.L.; Zhu, J.Z. The Finite Element Method: Its Basis and Fundamentals; Butterworth-Heinemann: London, UK, 2005; Volume 1. [Google Scholar]
  46. Ostachowicz, W.M.; Güemes, A. New Trends in Structural Health Monitoring; Springer: Berlin, Germany, 2013. [Google Scholar]
  47. Ván, P.; Berezovski, A.; Fülöp, T.; Gróf, G.Y.; Kovács, R.; Lovas, Á.; Verhás, J. Guyer-Krumhansl–type heat conduction at room temperature. EPL 2017, 118, 50005. [Google Scholar]
  48. Both, S.; Czél, B.; Fülöp, T.; Gróf, G.Y.; Gyenis, Á.; Kovács, R.; Ván, P.; Verhás, J. Deviation from the Fourier law in room-temperature heat pulse experiments. J. Non-Equilib. Thermodyn. 2016, 41, 41–48. [Google Scholar] [CrossRef]
  49. Kovács, R.; Ván, P. Models of ballistic propagation of heat at low temperatures. Int. J. Thermophys. 2016, 37, 95. [Google Scholar] [CrossRef]
  50. Kovács, R.; Ván, P. Generalized heat conduction in heat pulse experiments. Int. J. Heat Mass Transf. 2015, 83, 613–620. [Google Scholar] [CrossRef]
  51. Hofer, M.; Finger, N.; Kovacs, G.; Schoberl, J.; Zaglmayr, S.; Langer, U.; Lerch, R. Finite-element simulation of wave propagation in periodic piezoelectric SAW structures. IEEE Trans. Ultrason. Ferroelectrics Freq. Control 2006, 53, 1192–1201. [Google Scholar] [CrossRef]
  52. Rieth, A.; Kovács, R.; Fülöp, T. Implicit numerical schemes for generalized heat conduction equations. Int. J. Heat Mass Transf. 2018, 126, 1177–1182. [Google Scholar] [CrossRef]
  53. Kudela, P.; Krawczuk, M.; Ostachowicz, W. Wave propagation modelling in 1D structures using spectral finite elements. J. Sound Vib. 2007, 300, 88–100. [Google Scholar] [CrossRef]
Figure 1. Some typical two-dimensional B-spline wavelets on interval BSWI43 scaling functions: (a) corner function; (b) boundary function; (c) inner function; (d) partly inner function.
Figure 1. Some typical two-dimensional B-spline wavelets on interval BSWI43 scaling functions: (a) corner function; (b) boundary function; (c) inner function; (d) partly inner function.
Materials 12 01337 g001
Figure 2. Thermal wave propagation in one-dimensional media: (a) λ T = 0 ; (b) λ T = 0.1 ; (c) λ T = 0.5 ; (d) λ T = 1.5 .
Figure 2. Thermal wave propagation in one-dimensional media: (a) λ T = 0 ; (b) λ T = 0.1 ; (c) λ T = 0.5 ; (d) λ T = 1.5 .
Materials 12 01337 g002
Figure 3. Performance and convergency of finite element method (FEM) on one-dimensional issue (a) X = 0; (b) X = 0.6283; (c) X = 1.257.
Figure 3. Performance and convergency of finite element method (FEM) on one-dimensional issue (a) X = 0; (b) X = 0.6283; (c) X = 1.257.
Materials 12 01337 g003
Figure 4. Comparisons on accuracy and efficiency among FEM, spectral finite element method (SFEM), and wavelet finite element (WFEM): (a) X = 0; (b) X = 0.6283; (c) X = 1.257.
Figure 4. Comparisons on accuracy and efficiency among FEM, spectral finite element method (SFEM), and wavelet finite element (WFEM): (a) X = 0; (b) X = 0.6283; (c) X = 1.257.
Materials 12 01337 g004
Figure 5. WFEM performances for different λ T : (a) X = 0; (b) X = 0.6283; (c) X = 1.257.
Figure 5. WFEM performances for different λ T : (a) X = 0; (b) X = 0.6283; (c) X = 1.257.
Materials 12 01337 g005
Figure 6. Mesh grids used in calculations: (a) grid A; (b) grid B.
Figure 6. Mesh grids used in calculations: (a) grid A; (b) grid B.
Materials 12 01337 g006
Figure 7. Thermal behaviors (snapshots) obtained by mesh grid in Figure 6a with different λ T : (a) wavy behavior λ T = 0 ; (b) wavelike behavior λ T = 0.1 ; (c) diffusive behavior λ T = 0.5 ; (d) over diffusive behavior λ T = 1.5 .
Figure 7. Thermal behaviors (snapshots) obtained by mesh grid in Figure 6a with different λ T : (a) wavy behavior λ T = 0 ; (b) wavelike behavior λ T = 0.1 ; (c) diffusive behavior λ T = 0.5 ; (d) over diffusive behavior λ T = 1.5 .
Materials 12 01337 g007
Figure 8. Thermal behaviors (snapshots) obtained by mesh grid in Figure 6b with different λ T : (a) wavy behavior λ T = 0 ; (b) wavelike behavior λ T = 0.1 ; (c) diffusive behavior λ T = 0.5 ; (d) over diffusive behavior λ T = 1.5 .
Figure 8. Thermal behaviors (snapshots) obtained by mesh grid in Figure 6b with different λ T : (a) wavy behavior λ T = 0 ; (b) wavelike behavior λ T = 0.1 ; (c) diffusive behavior λ T = 0.5 ; (d) over diffusive behavior λ T = 1.5 .
Materials 12 01337 g008
Figure 9. WFEM performances for different λ T : (a) (X, Y) = (−1.508, 0); (b) (X, Y) = (−0.754, 0); (c) (X, Y) = (0, 0).
Figure 9. WFEM performances for different λ T : (a) (X, Y) = (−1.508, 0); (b) (X, Y) = (−0.754, 0); (c) (X, Y) = (0, 0).
Materials 12 01337 g009
Table 1. The main algorithm.
Table 1. The main algorithm.
Main Algorithm
#1 Loop over elements e
  Calculate the elemental matrices Ke, Me, and Ce
  Assemble matrices M, C, vector G
  Store every elemental stiffness matrix K
#1 End of loop over element e
Calculate the auxiliary vectors M 0 , M 1 , and M 2
Apply the initial condition
#2 Loop over time instants λ
  #21 Loop over elements e
Load the stiffness matrix Ke
Calculate f ^ e = K e Θ λ e on elemental level
Assemble vector F ^ by f ^ e
  #21 End of loop over elements e
Calculate effective vector R ˜ = G λ F ^ + M 1 Θ λ M 2 Θ λ Δ λ
Calculate Θ λ + Δ λ = ( M 0 ) 1 R ˜
#2 End of loop over time instants λ .

Share and Cite

MDPI and ACS Style

Yang, Z.-B.; Wang, Z.-K.; Tian, S.-H.; Chen, X.-F. Analysis and Modelling of Non-Fourier Heat Behavior Using the Wavelet Finite Element Method. Materials 2019, 12, 1337. https://doi.org/10.3390/ma12081337

AMA Style

Yang Z-B, Wang Z-K, Tian S-H, Chen X-F. Analysis and Modelling of Non-Fourier Heat Behavior Using the Wavelet Finite Element Method. Materials. 2019; 12(8):1337. https://doi.org/10.3390/ma12081337

Chicago/Turabian Style

Yang, Zhi-Bo, Zeng-Kun Wang, Shao-Hua Tian, and Xue-Feng Chen. 2019. "Analysis and Modelling of Non-Fourier Heat Behavior Using the Wavelet Finite Element Method" Materials 12, no. 8: 1337. https://doi.org/10.3390/ma12081337

APA Style

Yang, Z. -B., Wang, Z. -K., Tian, S. -H., & Chen, X. -F. (2019). Analysis and Modelling of Non-Fourier Heat Behavior Using the Wavelet Finite Element Method. Materials, 12(8), 1337. https://doi.org/10.3390/ma12081337

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop