Introduction to Numerical Approaches for Forward and Inverse Heat Transfer Problems

1. Introduction

Heat flow is described by the heat (or diffusion) partial differential equation, which takes the form ∂u∂t=k∂2u∂x2with u=uxtrepresenting the temperature at time tand at spatial location x(which may have one or more components). The constant kis the ratio of the thermal conductivity in the material over the specific heat times the density of the material. Applying the conservation law principle to the heat flow across the rod yields the parabolic partial differential equation. If a heat source is present or the material is exothermic, a non-homogeneous version of the problem: ∂u∂t−k∂2u∂x2=fxtwith fxt≠0applies. The initial boundary value problem comes from assigning relevant initial and boundary conditions. For example, for a perfectly insulated rod as in Figure 1, the conditions could be u0t=uLt=0. Given some initial profile of the temperature at time t=0, ux0=hx, the analytical solution in 1D or 2D with a regular geometrical domain is given by separation of variables. Assuming uxt=FxGt, one obtains a system of two ODEs and from that, the series solution uxt=∑n=1∞unxt=∑n=1∞Anexp−λn2tsinnπLxwith the coefficients Ansatisfying the initial condition ux0=∑n=1∞AnsinnπLx=fxand the eigenvalues λn=nkπL. The coefficients are evaluated via the sin series integration: An=2L∫0LfxsinkπxLdx. For non-intregrable cases or more complicated geometries, numerical techniques for approximation are used [1].

1.1 Derivation of homogeneous and non-homogeneous equations

We now discuss two simple scenarios which can be used to derive homogeneous and non-homogeneous forms of the heat Equation [2]. Heat is the flow of energy from a warmer to a cooler location. There are several types of heat transfer including conduction (flow of heat through stationary material), convection (flow of heat through fluids), and radiation (flow of heat through electromagnetic waves) [2]. In this chapter we are concerned mainly with heat conduction for which the differential heat model applies. We take Qto represent the heat transfer rate in units of Watts. It is the rate at which heat is flowing through one of the slabs in Figure 2. When heat flows into and out of a mass, the temperature of the mass changes in accordance with the model: Qin−Qout=mcdTdt, where Trepresents the temperature, mthe mass, and cis the specific heat value of the material. It should be noted that in heat problems, Tand uare ofter used interchangeably. We define q=Q/A(heat over area) as the heat flux, which by Fourier’s law is proportional to the temperature gradient: q=−k¯dTdx, where k¯is the thermal conductivity of the material [2]. Using Qin=Qxand Qout=Qx+δxalong with heat flux q=Q/A, we have Aqx−Aqx+δx=mc∂T/∂t. Applying Fourier’s law yields:

k¯A−∂T∂xx+∂T∂xx+δx=mc∂T/∂t

Next, using Taylor’s approximation ∂T∂xx+δx≈∂T∂xx+∂2T∂x2xδxand expressing m=ρAδxin terms of slab density, area, and width, we obtain:

k¯A∂2T∂x2=ρc∂T∂t⇒k∂2T∂x2=∂T∂t

For the exothermic slab case exhibited in Figure 2 on the right, we obtain an inhomogeneous equation. We assume the slab is filled with a material which emits heat at a rate q̇with units W/m3(energy per unit volume). In this case, the in and out energy balance gives Qx+q̇V−Qx+δx=mc∂T/∂t, which with q=Q/A,V=Aδx, Fourier’s law and Taylor’s expansion, gives k¯A∂2T∂x2δx+q̇Aδx=mc∂T∂t⇒k¯∂2T∂x2+q̇=ρc∂T∂t⇒k∂2T∂x2+q̇ρc=∂T∂t. In the case of steady state, ∂T/∂t=0and a second order differential equation model results.

2. Forward problem

The forward problem is concerned by the prediction of heat profile in future time from the initial condition. We first consider the heat equation initial boundary value problem in 1D:

with constant k∈ℝ. Via a transformation, this problem can be extended to nonzero boundary conditions. For (3), separation of variables gives the solution as the superposition of solutions of the form:

A quick calculation shows that ∂∂t−k∂2∂x2u˜xt=0. In order to form the well-known explicit scheme, we use the forward in time and central in space difference operators:

which results in the explicit finite difference scheme:

This can be re-written as ui,j+1=1−2αui,j+αui+1,j+ui−1,jwith α=klh2, with lthe time step (Δt) and hthe spatial step size (Δx). Von Neumann stability analysis yields the condition l≤h22k. For non-homogeneous cases the general model is ∂u∂t=k∂2u∂x2+fxtand the finite difference analysis proceeds in a similar way, with an explicit scheme of the form: ui,j+1=ui,j+kl/h2ui+1,j−2ui,j+ui−1,j+lfxitj.

If we use the backward time finite difference approximation for ∂u∂t, we get: ui,j+1−ui,jk=ui+1,j+1−2ui,j+1+ui−1,j+1h2, yielding the implicit scheme: −αui−1,j+1+1+2αui,j+1−αui+1,j+1=ui,jwhich can be solved for time j+1from the information at time jusing a linear system with a tridiagonal matrix with elements −α1+2α−α. The implicit scheme is numerically stable for any αbut is costlier than the explicit scheme. Another example, is the so called Crank–Nicolson scheme based on the time derivative approximation at the midpoint of the time stamp. It takes the form:

Again using α=klh2, we arrive at the scheme:

Note that this takes the form of the linear system Au=b, where Ais a tridiagonal matrix of terms −α21+α−α, as before, ucontaining the heat distribution at time j+1and bcontaining the heat distribution information at time j. A linear system solves then yields the new time information from previous data.

In 2D, the standard explicit scheme arises from the discretization:

Like in the 1D scheme, suitable finite difference approximations yield the implicit method.

A practically useful direction to consider is the use of non-standard finite difference (nsfd) schemes [5]. Based on (5), we define the discretized operators:

and proceed to evaluate d¯t−kld¯2xh2u˜xtfor u˜xt=e−π2ktsinπx. The standard fd scheme does not set this expression to zero, and hence this analytically valid solution does not satisfy the original partial differential equation. However, based on the calculations below, we can make a replacement for klh2, which would force the equation to zero.

Setting d¯t−klh2d¯2xu˜xt=0, we obtain:

Hence, the nsfd scheme follows if we make the following replacement in (6):

When h,l≪1, the nsfd substitution coincides with the standard finite difference scheme. We use the small ∣q∣Taylor series approximations eq≈1+qand sinq≈qto obtain:

By this it follows that for small h,l, the stability constraint klh2≤1is satisfied in the nsfd scheme when l≤h22just like in the original scheme. This also holds for large h, as long as his chosen so that ∣cosπh−1∣is not very small.

2.1 Numerical experiments

Below, we present some sample results for the case of constant k∈ℝin (3), where we found that the described nsfd scheme often outperforms the fd scheme at the same time step. In Figure 3, we illustrate the solutions with fd and nsfd methods with the initial condition, ux0=sinπx+xfor problem 1, with boundary conditions u0t=0and u1t=1and with the initial condition ux0=12sin9πx−7sin4πxand zero homogeneous boundary conditions for problem 2. For the second example, errors vs. step size are also compared. In general, we observe more accurate results with the use of the nsfd scheme.

An interesting case is that of heat transfer with an exothermic material, which can be modeled with a heat source distributed throughout [3]. Considering an infinitesimal segment as in Figure 1 yields the energy balance equation: q̇+αdx−q¯+dq̇dxdx=0⇒dq̇dx=α. It follows that the heat model takes the form 1α∂u∂t=∂2u∂x2+q̇k. The steady state behavior does not depend on time so that ∂u∂t=0and an ordinary differential equation model d2udx2+αk=0results. Upon a double integration, the general solution is ux=−α2kx2+Ax+B, which with the boundary conditions T0=TL=Twbecomes ux=−α2kx2+Lx+Tw, a parabola. We can simulate this steady state example with a finite difference scheme for the time derivative, yielding the temperature profile shown in Figure 4.

3. Inverse problem

The inverse problem consists of predicting the initial condition from the observation of the heat profile at some later time t. In general, the problem is ill-posed and regularization is necessary for its solution [4]. Let us look at a simple example in 1-D. Starting from the initial boundary value problem on x∈0L:

In this case, separation of variables can be used and after a well-known sequence of steps, the following series solution can be derived [3]:

It follows that ux0=∑n=1∞AnsinnπxLwhich is a sin series; that is, An=Fux0=2L∫0Lux0sinnπxL. It follows that the sin coefficients of the solution at t=Tare given by the sin coefficients at t=0multiplied by exp−knπL2T. The values of uat t=Tcan then be computed via the inverse sin transform applied to the product. In a similar way, it follows that from the transform of the solution at t=T, FuxT, we can in principle go backwards, to obtain ux0=F−1FuxT.∗expknπL2Twith the .∗representing component-wise operation [3]. However, with the negative sign inside the exponential gone, there is potential for blow up of any large magnitude terms. The issue can be addressed by performing filtering on the large coefficients, setting the largest magnitude portion to zero or some maximum value. The magnitude cut off value can be selected based on quartile statistics of the nonzero absolute values.

A more general approach for inversion is to utilize the inverse problem formulations of the matrix heat propagation formulation. For example, for the standard implicit or Crank-Nicholson scheme, we have the system Au=b, where Ais a tridiagonal matrix of terms −α21+α−α. This means that u1=A−1bu0, u2=A−1bu1=A−2bu0⇒bu0=A2u2. This then implies that the initial condition at time t=0, u0can be recovered from the application of the inverse matrix raised to a power ntcorresponding to the time point. As another example, if we discretize ut=kuxxas:

with M=kh2tridiag−12−1, for j=0,1,…,nt−1. Using ux0=U0, we obtain the relations:

This means that U0can in theory be obtained right from I+ΔtMntunt. The matrix needs not be formed explicitly, as the finite differencing scheme can be repeatedly applied to unt=uT. However, as shown in Figure 5, raising the typical differencing matrix to a power results in non-linear singular value decay. Hence, the application of the inverse of matrix I+ΔtM−ntto the approximant for b=uT, the heat profile at time T, would result in the blow-up of any numerical errors due to the ill-conditioning. We must instead consider the inverse formulation LU0=bwith L=L¯−ntand utilize a regularization scheme.

This can be accomplished by replacing the naive solution L−1bby the regularized formulation minu∥Lu−b∥22+λ∥u∥22, which for large enough λ>0, alleviates the problem due to the fast decay of singular values. This well-known Tikhonov regularization formulation essentially acts to remove the influence of small singular values on the solution, which would otherwise blow up in value upon inversion. The regularization parameter λshould be chosen with care to provide regularization but not to move too far away from the naive inverse operator. For this reason, λis usually selected based on an L-curve criterion [6], whereby the parameter is initially set at some fraction of ∥LTb∥∞and then iteratively lowered to a lower fraction (with the values either linearly or logarithmically spaced). At each λvalue, the solution to LTL+λIuλ=Ltbis obtained (e.g. with a Conjugate Gradient scheme) using the solution at the previously used λas the initial guess. Typically, the progression to lower λstops once the value of ∥Luλ−b∥2/Nreaches a plateau. It is possible to employ a predictor-corrector scheme to iteratively improve the solution. In this approach, the approximation to U0can be obtained as above, then the forward operator can be applied to propagate this approximated initial heat profile at t=0to time t=Tand compare with observations [4]. We can compare the residual at time Tand attempt some corrections to U0based on this residual such that the propagated solution at t=Tbest matches observations. In the case of only a finite subset of observations, interpolation techniques can be utilized to complete the profile at t=Tacross all values of x[7] prior to regularization to obtain an approximate smooth and continuous uTheat profile.

3.1 Numerical experiments

As an example, consider the 1D cosine based initial heat profile with k=1and zero boundary conditions. The profile and subsequent integration results are shown in Figure 6. It is clear that the temperature in the rod will approach zero at all locations with increasing time. The inverse problem consists of approximately recovering the initial profile ux0=0.5cos2π/65x+1), with support on N/43N/4. For this, we take the sin transform of the solution at t=Tand multiply by expknπL2T. In Octave, this can be accomplished with the code:

fuT=fft (u (:, n)); % temperature profile at time T

R=exp (k.*k ′*T);

filter (R);

init_cond_approx=real (ifft (fuT.*R));

The term Rmay contain large entries, which will blow up when the inverse fft is performed. For this reason, the filtering of small entries is necessary. This can be accomplished simply by replacing all large entries of Rwith either zero or some max value; e.g. taking the cutoff to be all values above 1e6and replacing them either with 0 or 1e2. The cutoff can be estimated with standard statistical outlier techniques, by taking the quartiles Q1,Q2,Q3of the ∣R∣values, computing IQR=Q3−Q1and filtering all entries of Rabove Q3+αIQRwith α≥2by setting them to the Q2value. The process is illustrated in Figure 7, for the same example as above with filtering for absolute value of entries above 1e1,1e3,1e5. It can be observed that filtering only the entries of Rabove 1e5in absolute value leads to a good approximation of the initial condition. Tikhonov regularization can be utilized when it is possible to build up the finite difference propagation operator (L), either as an explicit matrix, or as a function. An example is shown to the right of Figure 7, with an initial triangular temperature profile that is approximately recovered by a simple continuation scheme. In the same Figure, we also show the profiles of uTboth as observed and as propagated from the approximated initial values. The smoothing properties of the heat kernel [3] give rise to the observed transformation between t=0and t=T.

4. Conclusions

This article discusses numerical approaches for the heat equation, both for the forward and inverse problems with examples for homogeneous and non-homogeneous cases. For the forward heat propagation modeling problem, we have presented a simple non-standard finite difference formulation, which satisfies the differential model and appears to result in improved performance over the standard finite difference scheme with similar magnitude time and spatial step sizes. For the inverse problem, we have described approaches based on Fourier filtering and Tikhonov regularization.