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

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.


Introduction
The classical Fourier heat conduction law follows a property, which is not physical, that is, an infinite thermal wave propagation speed [1][2][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]: 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][8][9][10]. Then the required moments are found from the coefficient of Taylor's series [9][10][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 2 The Scientific World Journal 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).

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 where 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 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 ..
The Laplace transform of (5) is shown below: The system solution ( ) can be approximated by using polynomial equation in s-domain, as given by the following equation: The moments, , are the coefficients of Taylor series expansion about = 0 (Maclaurin series).

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.
The Scientific World Journal 3 3.1. Zero State Response. In ZSR, the initial condition is assumed to be zero, (0) = 0 and (0) = 0: By equating the same powers of , the moments are generated from the equation below: where is the order of Pade approximation.

Zero Input Response.
In ZIR, the input force is assumed to be zero, = 0: By equating the same powers of , the moments are generated from the equation below: 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 ‖ − ‖ 2 becomes smaller for the proper choice of = ( * ) −1 * [16]. The singularity condition of * is removed by adding ℎ 2 term where, after modification, it becomes 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 ℎ = ( * + ℎ 2 ) −1 . 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 = ( * 1 , * 2 , . . . , * ) and the TWCAWE moments subspace is denoted as follows: * = ( 1 , 2 , . . . , ), where the relation between ICAWE and TWCAWE moments is 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): The nodal moment [ ] can be extracted from the global moment matrix [ ] for any arbitrary node as follows: 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 Poles and residues can be found by solving ] . 4 The Scientific World Journal

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 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)

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

Transient Response. The final transient response for any arbitrary node is
Instantaneous temperature rise

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 The Scientific World Journal . 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. 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.
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.
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.

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.