Fast Transient Thermal Analysis of Non-Fourier Heat Conduction Using Tikhonov Well-Conditioned Asymptotic Waveform Evaluation

Abstract

Non-Fourier heat conduction model with dual phase lag wave-diffusion model was analyzed by using well-conditioned asymptotic wave evaluation (WCAWE) and finite element method (FEM). The non-Fourier heat conduction has been investigated where the maximum likelihood (ML) and Tikhonov regularization technique were used successfully to predict the accurate and stable temperature responses without the loss of initial nonlinear/high frequency response. To reduce the increased computational time by Tikhonov WCAWE using ML (TWCAWE-ML), another well-conditioned scheme, called mass effect (ME) T-WCAWE, is introduced. TWCAWE with ME (TWCAWE-ME) showed more stable and accurate temperature spectrum in comparison to asymptotic wave evaluation (AWE) and also partial Pade AWE without sacrificing the computational time. However, the TWCAWE-ML remains as the most stable and hence accurate model to analyze the fast transient thermal analysis of non-Fourier heat conduction model.

1. Introduction

The classical Fourier heat conduction law follows a property, which is not physical, that is, an infinite thermal wave propagation speed [1–3]. Fourier’s law is quite accurate for most common engineering problems. In some cases, the Fourier law cannot explain, such as near the absolute zero temperature, thermal gradient, which is extreme [3]. Non-Fourier conduction accounts for two-phase lag due to thermal wave propagation speed and another one due to relaxation time of electron. D. Y. Tzou proposed a non-Fourier model based on two-phase lag as follows [1, 3, 4]: 𝑄𝑟,𝑡+𝜒_𝑞=−𝐾_𝑐∇𝜃𝑟,𝑡+𝜒_𝑇, (1) where 𝜒_𝑇 stands for finite relaxation time for electron phonon and 𝜒_𝑞 represents finite thermal wave speed [5].

Conventional solver based on iterative, as fourth-order Runge-Kutta technique, is computationally expensive, though its result is accurate. In contrast, asymptotic waveform evaluation (AWE) technique is much faster solver [6]. Taylor’s series was used to expand the AWE transfer function [7–10]. Then the required moments are found from the coefficient of Taylor’s series [9–11]. AWE technique is able to solve up to three-dimensional models [12]. However, the limitation of AWE model is that this model cannot predict the temperature responses accurately as it is ill-conditioned moment matching [6, 13]. Then Loh et al. [6] developed technique, which is based on the AWE to remove the instability of AWE, called partial Pade AWE. In partial Pade AWE, selected poles and residues from any heat imposed boundary are used to calculate the temperature responses of the whole system. However, this technique is unable to predict the actual temperature. This prompted the present study to use another technique, which is based on well-conditioned asymptotic waveform evaluation (WCAWE).

WCAWE algorithm was introduced by Slone et al. [13, 14] to solve frequency domain electromagnetic problem based on the FEM. In this model, a correction term was introduced during the poles and residues calculation. Maximum likelihood (ML) was used successfully to calculate poles and residue to bring down the instability of AWE and also partial Pade AWE. Later, Liu et al. solved Fourier heat conduction by using the WCAWE technique [15]. The model proposed in the present work is based on the WCAWE algorithm to analyze the non-Fourier thermal problem. To avoid the singularity problem, our proposed model is embedded with Tikhonov regularization technique [16] to further enhance stable responses. Hence, the present model is recognized as TWCAWE. Later, we developed another well-conditioned scheme, called TWCAWE using the mass effect (TWCAWE-ME) that reduces the computational time compared to TWCAWE using the ML scheme (TWCAWE-ML).

2. Non-Fourier Model

The non-Fourier model with two-phase lag is shown in (2) in the case of fast thermal conduction. This model is converted into a normalized hyperbolic equation given by 𝜕²𝜑𝜕𝜍²+𝜕²𝜑𝜕𝜓²+𝑍_𝑇𝜕³𝜑𝜕𝜏𝜕𝜍²+𝑍_𝑇𝜕³𝜑𝜕𝜏𝜕𝜓²=𝜕𝜑𝜕𝜏+𝑍_𝑞𝜕²𝜑𝜕𝜏², (2) where 𝜑=𝜃−𝜃₀𝜃_𝑤−𝜃₀𝑡,𝜏=𝑙²/ℎ⁻𝑥,𝜍=𝑙,𝑦𝜓=ℎ,𝑍_𝑇=𝜒_𝑇𝑙²/ℎ⁻,𝑍_𝑞=𝜒_𝑞𝑙²/ℎ⁻, (3) where ℎ⁻ is the thermal diffusivity.

FEM meshing is performed to generate the rectangular element. The nodes of the element are denoted by 𝑗, 𝑘, 𝑚, and 𝑜. Applying FEM based on the Galerkin weighted residual technique of (2), we obtain _Δ[𝑁]^𝑇𝜕^2𝜑𝜕𝜍²+𝜕²𝜑𝜕𝜓²+𝑍_𝑇𝜕³𝜑𝜕𝜏𝜕𝜍²+𝑍_𝑇𝜕³𝜑𝜕𝜏𝜕𝜓²−𝜕𝜑𝜕𝜏−𝑍_𝑞𝜕²𝜑𝜕𝜏²𝜕𝜍𝜕𝜓=0. (4)

The above non-Fourier heat conduction equation shown in (4) is converted into a linear equation based on the Galerkin weighted residual technique, as shown in 𝐷𝜃^..+𝑉𝜃̇+𝐾𝜃=𝐹(𝑡), (5) where mass matrix 𝐷∈𝑅^𝑁×𝑁, damping matrix 𝑉∈𝑅^𝑁×𝑁, stiffness matrix 𝐾∈𝑅^𝑁×𝑁, and 𝐹(𝑡)∈𝑅^𝑁. The elemental matrices obtained from (4) are as follows: 𝐷=𝐴𝑍_𝑞,𝜍364212242112422124𝐾=+𝜓6∗𝜓21−1−212−2−1−1−221−2−112,𝐴6∗𝜍21−1−212−2−1−1−221−2−112𝑉=+𝑍364212242112422124^𝑇𝜍+𝑍6∗𝜓21−1−212−2−1−1−221−2−112_𝑇𝜓.6∗𝜍21−1−212−2−1−1−221−2−112 (6) The Laplace transform of (5) is shown below: 𝐷𝑠²𝜃(𝑠)−𝑠𝜃^′(0)−𝜃^′(0)((𝑠)(0(𝑠)𝑠+𝑉𝑠𝜃−𝜃))+𝐾𝜃=𝐹,₂𝐷𝜃(𝑠)+𝑠𝑉𝜃(𝑠)−𝐷𝜃^′(0(𝑠)(0))+𝐾𝜃−𝐷𝜃−𝑉𝜃^′(0)=𝐹. (7) The system solution 𝑇(𝑠) can be approximated by using polynomial equation in s-domain, as given by the following equation: 𝜃(𝑠)=^∞_𝑛=0𝑀_𝑛𝑠^𝑛. (8) The moments, 𝑀_𝑛, are the coefficients of Taylor series expansion about 𝑠=0 (Maclaurin series).

3. TWCAWE Algorithm

The TWCAWE algorithm is initiated by applying the Laplace transform on (5) where the moments obtained are used to estimate the zero state response (ZSR) and zero input response (ZIR). The moment computation for ZSR and ZSR is shown in (10) and (12), respectively.

3.1. Zero State Response

In ZSR, the initial condition is assumed to be zero, 𝜃(0)=0 and 𝜃^′(0)=0: 𝐷𝑠²𝜃(𝑠)−𝑠𝜃^′(0)−𝜃^′(0)+𝑉(𝑠𝜃(𝑠)−𝜃(0))+𝐾𝜃(𝑠)=𝐹,𝐷𝑠²+𝑉𝑠+𝐾𝜃(𝑠)=𝐹,𝐷𝑠²+𝑉𝑠+𝐾𝑀₀+𝑀₁𝑠+𝑀₂𝑠²+⋯+𝑀_𝑛𝑠^𝑛=𝐹. (9) By equating the same powers of 𝑠, the moments are generated from the equation below: 𝐾𝑀₀=𝐹𝐾𝑀₁=−𝑉𝑀₀𝐾𝑀_𝑛=−𝐷𝑀_𝑛−2+𝑉𝑀_𝑛−1,for𝑛=2,3,…,2𝑞−1, (10) where 𝑞 is the order of Pade approximation.

3.2. Zero Input Response

In ZIR, the input force is assumed to be zero, 𝐹=0: 𝐷𝑠²𝜃(𝑠)−𝑠𝜃^′(0)−𝜃^′(0)((𝑠)(0(𝑠)𝑠+𝑉𝑠𝜃−𝜃))+𝐾𝜃=0,₂𝐷+𝑠𝑉𝜃(𝑠)−𝐷𝜃^′(0)+𝐾𝜃(𝑠)=𝐷𝜃^′(0)+𝑉𝜃(0)𝑠²𝐷+𝑠𝑉𝜃(𝑠)−𝐷𝜃^′(0)×𝑀+𝐾₀+𝑀₁𝑠+𝑀₂𝑠²+⋯+𝑀_𝑛𝑠^𝑛=𝐷𝜃^′(0)+𝑉𝜃(0). (11) By equating the same powers of 𝑠, the moments are generated from the equation below: 𝐾𝑀₀=𝑉𝜃(0)+𝐷𝜃̇(0),𝐾𝑀₁=𝐷𝜃(0)−𝑉𝑀₀,𝐾𝑀_𝑛=−𝐷𝑀_𝑛−2+𝑉𝑀_𝑛−1,for𝑛=2,3,…,2𝑞−1. (12)

The moment matching process in AWE is inherently ill conditioned. To overcome this instability, WCAWE was proposed [13]. In our proposed model, to reduce the instability, we removed the singularity problem by changing the stiffness matrix. Consider a well-conditioned approximation problem 𝐾𝑥≈𝑦; the residual ‖𝐾𝑥−𝑦‖₂ becomes smaller for the proper choice of 𝑥=(𝐾∗𝐾)⁻¹𝐾∗𝑦 [16]. The singularity condition of 𝐾∗𝐾 is removed by adding ℎ²_𝑐𝐼 term where, after modification, it becomes 𝑥=𝐾∗𝐾+ℎ²_𝑐𝐼⁻¹∗𝑦, (13) where ℎ_𝑐 is the regulation parameter, which depends on the order of the equation, whereas 𝐼 is the identical matrix. The approximate inverse family is defined by 𝐶_ℎ=(𝐾∗𝐾+ℎ²_𝑐𝐼)⁻¹𝐾. As given in (10) and (12), moment calculation involves inverse of stiffness matrix 𝐾. The resulting moments, especially higher order moments, are product of lower order moments and, hence, the singularity problem can accumulate the compounded error. Therefore, in TWCAWE, computing moments by the inverse of 𝐾 matrix is now replaced by inverse of 𝐶_ℎ. Subsequently, the ill-conditioned moments are converted to well-conditioned moments. The new moments, called TWCAWE moments, are obtained from ICAWE moments. The ICAWE moment subspace is 𝑀_𝑛=(𝑢^∗₁,𝑢^∗₂,…,𝑢^∗_𝑛) and the TWCAWE moments subspace is denoted as follows: 𝑀^∗_𝑛=(𝑢₁,𝑢₂,…,𝑢_𝑛), where the relation between ICAWE and TWCAWE moments is 𝑢_𝑛=𝑢^∗_𝑛/𝜂_𝑛 where 𝜂_𝑛=‖𝑢^∗_𝑛‖. Modified Gram-Schmidt technique is required to orthonormalize 𝑀^∗_𝑛 that gives a more accurate solution [11]. In TWCAWE algorithm, 𝑍 can be computed as shown in (16): Proj_𝑀∗𝑎𝑀^∗_𝑛=𝑀^∗_𝑛⋅𝑀^∗_𝑎𝑀^∗_𝑎,for𝑎=1,2,…,(𝑛−1),𝑢_𝑛=𝑀^∗_𝑛−Proj_𝑀∗𝑎𝑀^∗_𝑛,𝑀_𝑛=𝑍⁻¹𝑀^∗_𝑛. (14)(15)(16)

The nodal moment [𝑚] can be extracted from the global moment matrix [𝑀] for any arbitrary node 𝑖 as follows: 𝑚_𝑛𝑖=𝑀^∗_𝑛𝑖. (17) The transient response for any arbitrary node 𝑖 can be approximated by using Pade approximation and then further simplified to partial fractions [17], as shown in 𝜃_𝑖(𝑠)=𝑚₀+𝑚₁𝑠+𝑚₂𝑠²+⋯+𝑚_𝑛𝑠^𝑛=𝑑₀+𝑑₁𝑠+⋯+𝑑_𝑠−1𝑠^𝑠−11+𝑐₁𝑠+⋯+𝑐_𝑠𝑠^𝑠=𝑅₁𝑠−𝑃₁+𝑅₂𝑠−𝑃₂𝑅+⋯+^𝑠𝑠−𝑃_𝑠. (18)

Poles and residues can be found by solving 𝑚₀𝑚₁𝑚₂⋯𝑚_𝑞−1𝑚₁𝑚₂𝑚₃⋯𝑚_𝑞𝑚⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅_𝑞−1⋯⋯⋯𝑚_2𝑞−2𝑐_𝑠𝑐_𝑠−1⋅⋅𝑐₁𝑚=−^𝑠𝑚_𝑠−1⋅⋅𝑚_1𝑛𝑖=1,𝑐_𝑖𝑃^𝑖𝑃+1=0,₁₋₁𝑃₂⁻¹𝑃₃⁻¹⋯𝑃_𝑠⁻¹𝑃₁⁻²𝑃₂⁻²𝑃₃⁻²⋯𝑃_𝑠⁻²𝑃⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅_1−𝑠𝑃₂^−𝑠𝑃₃^−𝑠⋯𝑃_𝑠^−𝑠𝑅₁𝑅₂⋅⋅𝑅_𝑠𝑚=−⁰𝑚₁⋅⋅𝑚_𝑞−1. (19)

4. T-WCAWE Using Maximum Likelihood Scheme

In T-WCAWE model, a correction term 𝐽_𝑢 is introduced to obtain more stable responses. This correction term 𝐽_𝑢 is calculated from 𝑍 matrix. The 𝑍 matrix is computed from ill-conditioned moments and well-conditioned moments, as shown in (16). This 𝑍 matrix is nonsingular [13, 14]. The correction term 𝐽_𝑢 is obtained by using maximum likelihood (ML) from 𝐽_𝑢(𝑖)=^𝑟_𝑖=1𝑍⁻¹_{[𝑖∶𝑖+2𝑞+𝑤−1,𝑖∶𝑖+𝑤−1]}, (20) where 𝑟 is the number of nodes, 𝑤 is the order of the equation, and 𝑞 is the order of the product. Now, the new T-WCAWE poles and residues are calculated from (21) and (22), respectively: 𝑃^∗(𝑖)=𝐽_𝑢(𝑖)𝑃(𝑖),𝑅^∗(𝑖)=𝐽_𝑢(𝑖)𝑅(𝑖). (21)(22)

5. TWCAWE with Mass Effect Scheme

In partitioning of FEM, it was reported that the resulting interpolation function can lead to singular or ill-conditioned. This singularity occurs because some interpolations are linearly dependent and ill-condition occurs because the functions are very close to each other. To minimize this singularity and ill-condition, the mass matrix is changed by introducing a factor “𝜍” [18]. For small value of 𝜍, TWCAWE-ME scheme reduces the effect of mass, as given in 𝐷_𝑠=(1−𝜍)𝐷_𝑐+𝜍𝐷_𝐿, (23) where 0≤𝜍≤1.

Here, 𝐷_𝑐=𝜅Δ𝑥²,𝐷364212242112422124_𝐿=𝜅Δ𝑥²41111, (24) where 𝜅 is the grid aspect ratio and Δ𝑥 is the length of each element in 𝑥 axis. 𝐷_𝑐 and 𝐷_𝐿 are consistent and row sum-lumped mass matrices, respectively. The lumped mass matrix 𝐷_𝐿 is calculated using the total mass and allocating the lumped masses in the ratio of the diagonal element of the consistent mass matrix [19]. In this scheme, 𝐷 from (5) is replaced by new calculated 𝐷_𝑠.

5.1. Transient Response

The final transient response for any arbitrary node 𝑖 is 𝑇_𝑖(𝑡)=zsr(𝑡)+zir(𝑡), (25) where zsr(𝑡)=^𝑞_𝑖=1𝑅^∗_𝑖𝑃^∗_𝑖𝑒^{(real(𝑃∗𝑖𝑡)}×cosimag𝑃^∗_𝑖𝑡+sinimag𝑃^∗_𝑖𝑡,(𝑡)=−1zir^𝑞_𝑖=1𝑃^∗_𝑖𝑒^{(real(𝑃∗𝑖𝑡)}×cosimag𝑃^∗_𝑖𝑡+sinimag𝑃^∗_𝑖𝑡. (26)

6. Results Analysis

In the present study, the non-Fourier heat conduction is analyzed by considering two-dimensional rectangular configurations using FEM and Tikhonov based WCAWE algorithm. The simulation work was carried out on a rectangular two-dimensional slab meshed to generate the rectangular elements, as shown in Figure 1 where the instantaneous boundary condition was applied at the left boundary of the slab. Boundary condition employed for this simulation work is given as 𝜑=𝜃−𝜃₀𝜃_𝑤−𝜃₀, (27) where 𝜃_𝑤=rise of temperature at the left boundary edge from the initial value 𝜃_𝑜. The applied normalized boundary conditions are stated as follows: 𝜑0,𝜍,𝜓=1,𝜕𝜑𝜕𝜉1,𝜍,𝜓=0,𝜕𝜑𝜕𝜓𝜍,1,𝜓=0,𝜕𝜑𝜕𝜓𝜍,0,𝜓=0. (28)

Figure 1 

Rectangular two-dimensional slab.

Figures 2, 3, and 4 show the temperature distribution for the above rectangular slab considering different values of 𝑍_𝑇 as 0.5, 0.05, and 0.0001, respectively. The temperature spectrum shown in Figures 2–4 is plotted by considering the node situated in the middle of the slab; all the responses are calculated by TWCAWE-ML scheme. It is clear from the figure that all temperature responses are continuous and accurate for different values of 𝑍_𝑇. In this scheme, the individual poles and residues are used to calculate the individual temperature responses. Hence, this scheme accurately predicts the temperature distribution across the system. The initial high frequency and delay due to relaxation time of electron are clear in Figure 4 for 𝑍_𝑇=0.0001.

Figure 2 

Normalized temperature response spectrum along the center of the slab in the case of 𝑍𝑇=0.5.

Figure 3 

Normalized temperature response spectrum along the center of the slab in the case of 𝑍𝑇=0.05.

Figure 4 

Normalized temperature response spectrum along the center of the slab for the case of 𝑍𝑇=0.0001.

Though the fourth-order Runge-Kutta technique is slow, the results obtained from this method are more accurate due to its stability. Therefore, our proposed model is compared with the results obtained by Runge-Kutta technique in Figures 5 and 6.

Figure 5 

Comparison with TWCAWE-ML, TWCAWE-ME, Runge-Kutta technique, and AWE using partial Pade in the case of 𝑍𝑇=0.05.

Figure 6 

Comparison of TWCAWE-ML, TWCAWE-ME, Runge-Kutta technique, and ICAWE for the case of 𝑍𝑇=0.0001.

To show the accuracy of our proposed model, we compared both schemes, TWCAWE-ML and TWCAWE-ME, against AWE using partial Pade and Runge-Kutta technique, as shown in Figure 5. The comparison has been presented with 𝑍_𝑇=0.05. The temperature behavior obtained from both schemes of TWCAWE is similar to Runge-Kutta technique. Partial Pade AWE, however, is unable to predict the temperature behavior as Runge-Kutta technique. In partial Pade AWE, arbitrarily chosen poles and residues are used to approximate the temperature behavior for the whole system. Hence, this method was unable to calculate the temperature behavior for other nodes. On the other hand, both schemes of TWCAWE model use its own poles and residue to calculate the temperature behavior.

The comparison between our proposed two schemes is shown in Figure 6 where 𝑍_𝑇=0.0001. TWCAWE-ME is able to calculate the delay and high frequency, but the temperature behavior deviates from Runge-Kutta results. On the other hand, the temperature behavior obtained from TWCAWE-ML is in close proximity to Runge-Kutta technique.

To bring out the effectiveness of employing Tikhonov technique, we have performed another comparison in Figure 7 in the case of 𝑍_𝑇=0.0001. This figure shows that the temperature behavior obtained from AWE technique largely deviates from Runge-Kutta behavior. In AWE model, unstable poles make moment matching process ill conditioned.

Figure 7 

Comparison of TWCAWE-ML, WCAWE-ML, and ICAWE for 𝑍𝑇=0.0001.

The result of the WCAWE-ML model is not similar to Runge-Kutta technique due to the singularity problem of the stiffness matrix as shown in Figure 7. The error is compounded during moment calculation and it becomes more evident for higher order moment calculation. This causes WCAWE-ML to be unstable, thus producing inaccurate initial high frequency responses, making it difficult to observe the delay. However, by embedding Tikhonov technique into existing WCAWE-ML, singularity problem was minimized. This is carried out by adding a new term based on Tikhonov regulation technique as given in (13). The new term reduces the symmetrical nature of stiffness matrix, thus minimizing the singularity problem with negligible effect on the temperature behavior. By observing temperature behavior on node 5 which is near the imposed boundary condition, Tikhonov based WCAWE-ML (TWCAWE-ML) is in close proximity to Runge-Kutta technique, as shown in Figure 7, compared to WCAWE-ML which is devoid of Tikhonov technique. The model proposed in the present work is able to reduce the effect of singularity as well as instability of the previous model.

Table 1 shows the computational time required by analyzing different methods in the present work. From the table, TWCAWE-ML is found to be almost two times faster than Runge-Kutta technique. The orthonormalization and maximum likelihood approximation process in T-WCAWE poses computational load that makes it slower compared to AWE or partial Pade AWE. As TWCAWE-ME does not require orthonormalization and maximum likelihood approximation, it performs as fast as AWE or AWE using partial Pade despite some sacrifice in temperature response quality in comparison to TWCAWE-ML. However, TWCAWE-ME showed more stable and hence accurate results compared to either AWE or AWE using partial Pade.

Table 1 

Required time for different methods.

7. Conclusion

In the present work, the heat conduction based on non-Fourier model has been solved using FEM. The phase lag responsible for finite relaxation time is varied to analyze the capability of our proposed model for predicting the heat conduction in two-dimensional model. From the results, both schemes discussed in this work are able to reduce the instability of AWE. From the comparison, we can conclude that T-WCAWE using ML results is in close proximity to Runge-Kutta technique compared to TWCAWE-ME. However, TWCAWE-ME scheme is 3.3 times faster than the conventional solver, namely, Runge-Kutta technique (RK4), and needs 1.7 times less computational time than T-WCAWE with ML scheme. Despite this, the TWCAWE-ML is more accurate for predicting the initial frequency as well as delay based on the close proximity to Runge-Kutta technique, hence emerging as the most stable and accurate technique to analyze the fast transient thermal analysis of non-Fourier heat conduction model when compared to AWE, partial Pade AWE, and TWCAWE-ME.

Nomenclature

𝑅∶	Residue
𝑀:	Moment
𝑃:	Pole
𝜑:	Dimensionless temperature
𝜏:	Dimensionless time
𝜍:	Distance along “𝑥”
𝜓:	Distance along “𝑦”
𝑁:	Shape function matrix
𝐴:	Element area
𝜃𝑜:	Initial temperature
𝜃𝑤:	Raised temperature
𝑙:	Length
ℎ:	Width
𝑍𝑇:	Phase lag for temperature gradient
𝑍𝑞:	Phase lag for heat flux.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This research work is supported by the University of Malaya High Impact Research (HIR) Grant (no. 51) sponsored by the Ministry of Higher Education (MOHE), Malaysia, as well as Science Fund (SF010-2013) sponsored by Ministry of Science, Technology and Innovation (MOSTI), Malaysia.