Open access peer-reviewed chapter

Introduction to Numerical Approaches for Forward and Inverse Heat Transfer Problems

By Sergey Voronin

Submitted: September 23rd 2019Reviewed: September 8th 2020Published: October 27th 2020

DOI: 10.5772/intechopen.93928

Downloaded: 96

Abstract

We describe the formulation of forward and inverse problems in heat transfer and illustrate several techniques applicable to their solution.

Keywords

  • heat transfer
  • finite differences
  • regularization
  • Fourier transforms
  • interpolation

1. Introduction

Heat flow is described by the heat (or diffusion) partial differential equation, which takes the form ut=k2ux2with 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: utk2ux2=fxtwith fxt0applies. 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=1unxt=n=1Anexpλn2tsinLxwith the coefficients Ansatisfying the initial condition ux0=n=1AnsinLx=fxand the eigenvalues λn=nkπL. The coefficients are evaluated via the sin series integration: An=2L0LfxsinkπxLdx. For non-intregrable cases or more complicated geometries, numerical techniques for approximation are used [1].

Figure 1.

Perfectly insulated rod used a model of heat transfer, an exothermic slab, and sample modeled heat profile in 2D.

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: QinQout=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 AqxAqx+δx=mcT/t. Applying Fourier’s law yields:

Figure 2.

Heat conducting slabs without and with internal heat generation.

k¯ATxx+Txx+δx=mcT/t

Next, using Taylor’s approximation Txx+δxTxx+2Tx2xδxand expressing m=ρAδxin terms of slab density, area, and width, we obtain:

k¯A2Tx2=ρcTtk2Tx2=Tt

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̇VQx+δx=mcT/t, which with q=Q/A,V=Aδx, Fourier’s law and Taylor’s expansion, gives k¯A2Tx2δx+q̇Aδx=mcTtk¯2Tx2+q̇=ρcTtk2Tx2+q̇ρc=Tt. In the case of steady state, T/t=0and a second order differential equation model results.

1.2 Numerical schemes

Discretization of the derivative operators yields the finite difference formulation. A number of different schemes are possible based on the difference schemes used to replace the derivatives in the equation resulting in either an explicit method (where the update for uat the next time step is given by a direct calculation from the previous one) or an implicit method (where the update for uis accomplished via a linear solve). For 2D heat transfer, the model becomes α22ux2+2uy2=utand again, the finite difference method can be used for discretization. The initial value problems are either forward or inverse, in the sense that they either predict the temperature profile at later time from the knowledge of the initial temperature distribution at time t=0or in the inverse case, predict the temperature distribution at time t=0from knowledge of the temperature at a future time. Sometimes, the initial boundary value problem may be prescribed extra conditions, such as Neumann or other derivate condition type and the inverse problem could instead seek the initial values for such expressions from future data [3]. In general, the inverse problem is more challening than the forward problem due to being ill-posed. The discretized operator matrix which can be assembled from the difference formulation is not well conditioned and as a result the problem is overly sensitive to slight changes in input values. Regularization techniques need to be used for the solution [4].

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:

ut=kuxx,0<x<1,E1
ux0=fx,0<x<1,E2
u0t=u1t=0,t>0E3

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:

u˜xt=eπ2ktsinπx

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

utfdxitj=uxitj+1uxitjlE4
2ux2fdxitj=uxi+1tj2uxitj+uxi1tjh2E5

which results in the explicit finite difference scheme:

ui,j+1=ui,j+klh2ui+1,j2ui,j+ui1,jE6

This can be re-written as ui,j+1=12αui,j+αui+1,j+ui1,jwith α=klh2, with lthe time step (Δt) and hthe spatial step size (Δx). Von Neumann stability analysis yields the condition lh22k. For non-homogeneous cases the general model is ut=k2ux2+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,j2ui,j+ui1,j+lfxitj.

If we use the backward time finite difference approximation for ut, we get: ui,j+1ui,jk=ui+1,j+12ui,j+1+ui1,j+1h2, yielding the implicit scheme: αui1,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:

ui,j+1ui,jl=kui+1,j+12ui,j+1+ui1,j+1+ui+1,j2ui,j+ui1,j2h2

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

αui+1,j+1+21+αui,j+1αui1,j+1=αui+1,j+21αui,j+αui1,j

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:

Ti,jn+1Ti,jnΔt=kTi,j+1n2Ti,jn+Ti,j1nΔx2+Ti+1,jn2Ti,jn+Ti1,jnΔy2

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:

d¯tfxt=fxt+lfxtE7
d¯x2fxt=fx+ht2fxt+fxhtE8

and proceed to evaluate d¯tkld¯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.

d¯tu˜xt=sinπxeπ2kt+leπ2kt=sinπxeπ2kteπ2kl1=u˜xteπ2kl1E9
d¯x2u˜xt=eπ2ktsinπx+h2sinπx+sinπxhE10
=eπ2kt12iex+hex+h2sinπx+12iexhexhE11
=eπ2kt12ieiπxeiπheiπxeiπh4isinπx+eiπxeiπheiπxeiπhE12
=eπ2kt12ieiπxeiπh+eiπheiπxeiπh+eiπh4isinπxE13
=eπ2kt12i2cosπheiπxeiπx4isinπxE14
=eπ2kt12i4icosπhsinπx4isinπxE15
=2eπ2ktsinπxcosπh1=2u˜xtcosπh1E16

Setting d¯tklh2d¯2xu˜xt=0, we obtain:

u˜xteπ2kl12klh2u˜xtcosπh1=0eπ2kl1=2klh2cosπh1

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

klh2=12eπ2kl1cosπh1

When h,l1, the nsfd substitution coincides with the standard finite difference scheme. We use the small qTaylor series approximations eq1+qand sinqqto obtain:

12eπ2kl1cosπh1=eπ2kl12sin2πh2121π2kl1π2h22=klh2

By this it follows that for small h,l, the stability constraint klh21is satisfied in the nsfd scheme when lh22just like in the original scheme. This also holds for large h, as long as his chosen so that cosπh1is not very small.

2.1 Numerical experiments

Below, we present some sample results for the case of constant kin (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πx7sin4π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.

Figure 3.

Comparison of fd and nsfd scheme solutions for the 1-D heat equation for problems 1 and 2 above.

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̇+αdxq¯+dq̇dxdx=0dq̇dx=α. It follows that the heat model takes the form 1αut=2ux2+q̇k. The steady state behavior does not depend on time so that ut=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.

Figure 4.

Temperature distribution inside a thin slab with an exothermic material with homogeneous boundary conditions as simulated by an fd scheme.

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

ut=k2ux2E17
ux0=fx,u0t=0=uLtE18

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

uxt=n=1AnsinnπxLexpkL2t

It follows that ux0=n=1AnsinnπxLwhich is a sin series; that is, An=Fux0=2L0Lux0sinnπxL. It follows that the sin coefficients of the solution at t=Tare given by the sin coefficients at t=0multiplied by expkL2T. 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=F1FuxT.expkL2Twith 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=A1bu0, u2=A1bu1=A2bu0bu0=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:

uttkh2ui+1j+1t2uij+1t+ui1j+1t=0=uij+1uijΔtMuj+1t,

with M=kh2tridiag121, for j=0,1,,nt1. Using ux0=U0, we obtain the relations:

u0=U0;uj+1=I+ΔtM1uj=I+ΔtMj+1U0L¯j+1U0=uj+1

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+ΔtMntto 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.

Figure 5.

Tridiagonal structure of differencing operator matrix and singular value decay of operator raised to a power.

This can be accomplished by replacing the naive solution L1bby the regularized formulation minuLub22+λu22, 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 LTband 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λb2/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 expkL2T. In Octave, this can be accomplished with the code:

Figure 6.

Homogeneous 1D profile and change of temperature with time.

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 Rvalues, computing IQR=Q3Q1and 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.

Figure 7.

Initial condition reconstruction with different filtering levels and illustration of inversion for initial condition using Tikhonov regularization with parameter tuned via the L-curve approach.

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.

© 2020 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Sergey Voronin (October 27th 2020). Introduction to Numerical Approaches for Forward and Inverse Heat Transfer Problems, Inverse Heat Conduction and Heat Exchangers, Suvanjan Bhattacharya, Mohammad Moghimi Ardekani, Ranjib Biswas and R. C. Mehta, IntechOpen, DOI: 10.5772/intechopen.93928. Available from:

chapter statistics

96total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Applications of Heat Transfer Enhancement Techniques: A State-of-the-Art Review

By Suvanjan Bhattacharyya, Devendra K. Vishwakarma, Sanghati Roy, Ranjib Biswas and Mohammad Moghimi Ardekani

Related Book

First chapter

Introduction to Infrared Spectroscopy

By Theophile Theophanides

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us