Numerical Solution of Hyperbolic Heat Conduction Equation

. Hyperbolic heat conduction problem is solved numerically. The explicit and implicit Euler schemes are constructed and investigated. It is shown that the implicit Euler scheme can be used to solve eﬃciently parabolic and hyperbolic heat conduction problems. This scheme is unconditionally stable for both problems. For many integration methods strong numerical oscillations are present, when the initial and boundary conditions are discontinuous for the hyperbolic problem. In order to regularize the implicit Euler scheme, a simple linear relation between time and space steps is proposed, which automatically introduces suﬃcient amount of numerical viscosity. Results of numerical experiments are presented.


Introduction
The process of heat transfer usually is described by parabolic models. In the case, when heat conductivity coefficients depend on the solution itself the structure and properties of the solution can change significantly, we note the similarity with diffusion and convection models, when convection is a dominated process in some subregions (for some applications and numerical analysis results see [2]).
If the elapsed time of heat transport process is very small the inertion starts to be important and the wave nature of the heat energy transport becomes a dominant process. For example, the physical phenomena involved in ultrashort laser pulse interactions with solids are very complex. It is shown in [12] that the classical parabolic two temperature model, which is based on the Fourier law, is not sufficiently accurate. Because the laser pulse duration is very short, the relaxation time of free electrons in metals and the relaxation time in phonon collisions should be taken into account (see also our paper [1] and references given in it).

R. Čiegis
The following model was used in [1] to describe the interaction of ultrashort (picosecond and femtosecond) laser impulses with different materials: where T E and T L are temperatures of electrons and lattice, respectively, C E , C L are heat capacities per unit volumes, k E , k L are thermal conductivities, γ denotes the electron-lattice coupling factor. Several strategies can be applied to construct numerical algorithms for finding solutions of the hyperbolic heat transport models. In the first approach differential equations are discretized in space applying the finite element, Galerkin or finite difference methods, then some general computational tool for numerical solution of systems of hyperbolic equations is applied, e.g. explicit Mac-Cormack's predictor-corrector or Lax-Wendroff methods. Such algorithms and results of computational experiments are presented in [6,7] for one temperature hyperbolic heat conduction problems. The efficiency and accuracy of such algorithms are investigated by solving test problems. It is noted that for some test problems the numerical solutions are not monotone and numerical oscillations still remained even for very large numbers of grid points (see, [6]). In [14], a spectral method is used for approximation of the diffusion operator, the obtained system of ordinary differential equations is solved by explicit Runge-Kutta methods of different order. It is shown that the spectral approximation also results in a solution in which strong numerical oscillations are present. A special post-processing technique is proposed to reduce these oscillations. We note that a similar strategy was used to solve parabolic problems in [11], where implicit Runge-Kutta schemes were used to integrate the obtained system of ODEs.
The second approach, is to use numerical algorithms targeted for solution of systems of hyperbolic PDEs. For the numerical solution of one axisymmetric, hyperbolic heat conduction 2D problem with complex nonlinear boundary conditions second order explicit Total Variation Diminishing (TVD) schemes were used in [15].
The third approach is based on the fact that the relaxation time τ E,L of the temperature is small and in many applications it can be equal to zero in some subdomains, changing the type of the equation from hyperbolic to parabolic one. For many numerical integrators of the hyperbolic systems, the assymptotics τ → 0 makes the algorithm ill-posed. In order to obtain the effective and homogeneous solvers, the integration in time must be done by an algorithm which is effective for both hyperbolic and parabolic problems. Development of such algorithm is the main goal of this paper. We note that in [10] a standard Galerkin method is used to approximate equations in space and a Crank-Nicolson method is applied for integration in time. No stability and convergence analysis is given here, the performance of the proposed algorithm is verified by solving a 1D test case.
The rest of the paper is organized as follows. In Section 2, the parabolic and hyperbolic heat conduction problems are formulated. Both formulations are based on the same energy balance equation and only different constitutive relations are used to define equations for the flux. Discrete approximations of parabolic problem are investigated in Section 3. It is shown that these schemes coincide with explicit and implicit Euler schemes in the case of one parabolic equation. In Section 4 these schemes are modified in order to approximate the system of hyperbolic heat conduction problem. The stability is investigated by using spectral and energy methods. In Section 5 it is proved that a solution of the explicit Euler scheme is strongly oscillating if the initial and boundary conditions of the hyperbolic problem are discontinuous. It is shown that the implicit Euler introduces numerical viscosity and a relation between space and time steps is proposed in order to regularize the discrete scheme. Results of numerical experiments are presented where discrete solution is converging to the physical solution of the hyperbolic heat conduction problem and no numerical oscillations are presented. Some final conclusions are given in Section 6.

Problem Formulation
Let us consider the balance of heat in an arbitrary 1D subset (a, b) and apply the divergence theorem. We obtain the energy transport equation where e(x, t) is the density of internal energy, 1 0 e dx is the total amount of heat, q = q(x, t) is the heat flux. The internal energy density e depends on the temperature e = e 0 + ρcT , ρ is the mass density, c is the specific heat capacity.
Parabolic problem. In the classical heat transfer model, according to the Fourier law the heat flux is proportional to the temperature gradient [8]: is called the heat conductivity. Substituting both constitutive relations into (2.1) we obtain the parabolic initial -boundary value problem

R. Čiegis
Hyperbolic problem. As the laser pulse duration is very short, the relaxation time of free electrons in metals and the relaxation time in phonon collisions should be taken into account. The following model of modified heat flux can be used: here τ is the relaxation time (τ ≪ 1). Systems of delay differential equations are important in many areas of science and particularly in the biological sciences (e.g., population dynamics and epidemiology).
Expanding q(x, t + τ ) in Taylor's series, we get the following equation: Substituting this constitutive relations into (2.1) we obtain the hyperbolic initial -boundary value problem

Discrete Algorithms for Parabolic Problem
Our main goal is to develop robust and efficient numerical algorithms, which can be used for solution of both types of heat conduction problems avoiding explicit adaptation of algorithms to specific type of the equation. Such a situation is not new in construction of numerical algorithms for solution of PDEs. Examples of problems, which are described by different types of PDEs in different regions are given by: • Parabolic + Elliptic problems in different subregions for porous media flows [4,5].
• Diffusion + Convection problems with different domination is subregions.
Here we should take into account the Lax theorem, which defines boundaries for any attempt to build fast and robust solvers [8,9]. In general, it states that the approximation and stability of the discrete method are necessary properties for the convergence of the discrete solution to the solution of a differential equation.

Approximation of the parabolic problem
Here we restrict our analysis to two most popular algorithms, used to solve parabolic problems, i.e., the explicit and implicit Euler methods. For simplicity of notations, let us assume that ρ = 1, c = 1.
In this paper we use staggered space grids and construct numerical algorithms in a form convenient to solve the parabolic system of heat conservation (2.1)-(2.2).
The explicit Euler method. Let us define in [0, 1] two staggered grids Substituting the discrete equation for the flux into the first equation of (3.1), we get the well known explicit Euler finite difference scheme for solution of a classical parabolic heat conduction problem In this paper we use the standard notations of finite differences: The convergence analysis of this finite-difference scheme is well-known (see, e.g. [8,13]). For completeness of the presentation, here we will give only basic results. For sufficiently smooth solution of the differential problem, the approximation error of the explicit Euler finite difference scheme (3.1) is given by ψ ≤ C(h t + h 2 x ). The stability of this discrete scheme can be proved by using different methods (e.g., the maximum principle, spectral analysis or energy estimates [8,13]). If discrete steps satisfy the inequality then the convergence estimate follows from the Lax Theorem: where U ∞ = max 0≤j≤J |U j |.

R. Čiegis
The implicit Euler scheme. Conditional stability of the explicit Euler scheme can be improved if we use the implicit Euler algorithm for discrete integration in time: Substituting the equation for the flux into the first equation, we get the well known implicit Euler finite difference scheme For sufficiently smooth solution of the differential problem the approximation error of the implicit Euler scheme is given by The unconditional stability of this discrete scheme again can be proved by using different methods (e.g., the maximum principle, spectral analysis or energy estimates [8,13]). Thus the unconditional convergence estimate follows from the Lax Theorem:

Results of computational experiments
In order to compare qualitative properties of solutions of parabolic and hyperbolic heat conduction models, here we present results for a simple heat conduction problem The Neumann type boundary condition is approximated by the finite difference equation: In Figure 1 the profiles of solutions are presented at different time moments It follows from the presented results, that the given solution decreases monotonically with respect to x coordinate and increases monotonically with respect to time.

Discrete Schemes for the Hyperbolic Problem
In this section we generalize explicit and implicit Euler schemes in order to apply them for numerical solution of hyperbolic heat conduction problem (2.4). Again, two staggered grids are used to approximate differential equations in space.

The explicit Euler method
The non-stationary flux equation is approximated by the explicit scheme, thus we obtain the explicit Euler scheme for a system of differential equations (2.4): Substituting the discrete equation for the flux into the first equation, we get the following explicit Euler finite difference scheme for the hyperbolic heat conduction equation: (4.1)

R. Čiegis
Next we investigate the approximation and stability properties of the obtained finite difference scheme. This analysis is done for sufficiently smooth solution of the differential problem (2.4). Using Taylor's series it is easy to prove that the approximation error is given by Stability analysis. First, we use the spectral analysis, when k = const. The characteristic equation of the finite difference equation (4.1) is given by where the eigenvalues of the difference operator (−∂ x ∂xU ) are denoted by λ l ≤ 4/h 2 x . The Hurwitz stability criterion yields that |g l | ≤ 1 if We note from this estimate, that in the hyperbolic case the stability requirements for the explicit Euler scheme are weaker than in the parabolic case if τ ≫ h 2 x . Next we investigate the stability of the explicit Euler scheme by using the energy method. Here general diffusion coefficient k = k(x) is considered (see also [1]). The following well-known stability result is valid for three-level finite difference schemes written in a canonical form [13]: where the notation y n 0 t = (y n+1 − y n−1 )/(2h t ) is used.

Theorem 1. Let operators R and A be symmetric and positively definite
and they do not depend on t. Then the following conditions are sufficient for the stability of (4.3) with respect to the initial conditions where the norm is defined by Using a simple equality U t = U0 t + h t 2 Ut t , we obtain the canonical three level form of finite-difference scheme (4.1): i.e. operators in the canonical form are given by In order to use Theorem 1 it is sufficient to prove the estimate R > 1 4 A. Let us assume that k(z) ≤ k 1 , then 0 < A < 4k 1 h 2 z I. Thus the stability condition is satisfied if estimate k 1 h 2 t > (τ + 0.5h t )h 2 z is valid, which coincides with the relation (4.2).

The implicit Euler method
Next we consider the implicit Euler scheme, when the flux equation is approximated by fully implicit formula: Substituting the discrete equation for the flux into the balance equation, we get the implicit Euler finite difference scheme, which approximates the hyperbolic heat conduction equation: For sufficiently smooth solution of the system of differential equations (2.4) approximation error of finite difference scheme (4.4) is given by Stability analysis. First, the spectral stability analysis is done. The characteristic equation of finite difference scheme (4.5) is given by: It is easy to see, that in this case the Hurwitz stability criterion is satisfied unconditionally.
Next we apply the energy method to investigate the stability of the implicit Euler finite difference scheme. Its canonical three level form is defined as: i.e. in the canonical form operators B, R are given by After simple computations we conclude that all estimates of Theorem 1 are satisfied unconditionally and hence the implicit Euler scheme (4.5) is unconditionally stable in the energy norm.

Analysis of the Solution of Hyperbolic Heat Conduction Problem
It is well known that only the principal part of the linear partial differential equation is important in the analysis of characteristics of the equation (and classification of PDE). Thus let us consider a model hyperbolic problem Its characteristic curves are defined by straight lines x ± t/ τ /k = const. Using these characteristics, the second-order equation can be decoupled into a pair of scalar transport equations [8] ∂V ± ∂t ± ∂V ± ∂x = 0.
This analysis enables us to predict the velocity of the travelling wave solution for more general hyperbolic heat conduction problems. We expect that at time t a front of the right-moving wave will reach point x f (t) = t √ We solve a hyperbolic heat conduction problem [14] T (x, 0) = 0, q(x, 0) = 0, 0 ≤ x ≤ 1.
In numerical experiments we take τ = 0.1 and solve this problem using explicit Euler scheme. In the first experiment the number of points in the space grid is J = 100 and time step selected according stability requirements. The obtained results are given in Figure 2a. It follows that the wave moves with a predicted velocity, but the solution changes non monotonically in x coordinate and strong oscillations are presented.

The Gibbs phenomenon
In order to investigate the nature of these oscillations we have increased the number of grid points till J = 400. The results of computations are presented in Figure 2b and strong oscillations are still noticeable. This numerical effect arises due to discontinuous initial-boundary conditions in (5.1) and is known as Gibbs phenomenon [14]. It describes the peculiar manner in which the Fourier series of a piecewise continuously differentiable periodic function f behaves at a jump discontinuity. The Gegenbauer reconstruction procedure recovers spectral accuracy up to the discontinuity points (see, e.g. [14]), but this post-processing procedure cannot be applied for general multidimensional solvers.
In order to prove that numerical oscillations occurred due to discontinuous initial -boundary conditions, let us consider a case of smooth initial-boundary conditions  In Figure 3 two numerical solutions are shown. They are obtained by using J = 200 and J = 2000 grid points. In the first case the time step h t is too large in order to capture the continuous nature of the solution and numerical oscillations are noticeable again, but for sufficiently large J (and respectively, small time step h t ) these oscillations disappear.

Computational experiments: implicit Euler scheme
In this section we present results of computational experiments, when the implicit Euler scheme is used to solve the hyperbolic heat conduction problem (5.1). First, we have fixed the number of space grid points J = 200 and solved this problem with different values of time step h t = 0.001 and h t = 0.0004. Numerical results are presented in Figure 4a. The solutions are monotone and the stiffness of the front is increased for a smaller time step. But if time step h t is reduced further, numerical oscillations of the discrete solution are appearing again (see, Figure 4b).  The regularization effect of the implicit Euler scheme is preserved only for sufficiently large values of time step. It is obvious, that for very small time steps h t both schemes, i.e. explicit and implicit Euler integration algorithms, produce very similar and oscillating numerical solutions. We propose to use discrete steps h x and h t connected by the relation h t = Ch x . Results of computational experiments are presented in Figure 5. It follows from these results that numerical oscillations of the discrete solution are not apparent for the regularized implicit Euler scheme.
It is well known that for simple transport equation the regularization of the implicit Euler scheme and up-wind approximation is achieved due to numerical diffusion of non-symmetrical (only first-order accurate) discretization (see the analysis of modified differential equations in [9,17]). Thus the implicit Euler scheme automatically defines the "vanishing viscosity" solution, which is a correct physical solution.
For hyperbolic heat conduction problem the implicit Euler finite-difference scheme (4.4) (and (4.5) is defined on the staggered grids. Therefore the transport operators are approximated by the symmetrical finite difference scheme which is the second-order accurate in space. Asymmetry (and numerical diffusion) is introduced by the implicit Euler approximation in time coordinate. Let us recall the canonical form of this scheme: where the regularization operator R = τ I + 0.5h t (I + h t A) h 2 t . If h t = Ch x , then in R the diffusion term is of the same order as the fast relaxation term and therefore numerical oscillations of the discrete solution are damped.

Conclusions
It is shown that the implicit Euler scheme can be used to solve parabolic and hyperbolic heat conduction problems. The scheme is unconditionally stable for both problems. In order to regularize this scheme, a simple linear relation between time and space steps is proposed, which automatically introduces numerical viscosity and thus numerical oscillations are not appearing. Thus on the basis of implicit Euler scheme robust solvers of more general heat conduction problems can be constructed.