**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 *<sup>∂</sup><sup>u</sup> <sup>∂</sup><sup>t</sup>* <sup>¼</sup> *<sup>k</sup> <sup>∂</sup>*2*<sup>u</sup> <sup>∂</sup>x*<sup>2</sup> with *u* ¼ *u x*ð Þ , *t* representing the temperature at time *t* and at spatial location *x* (which may have one or more components). The constant *k* is 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 <sup>∂</sup><sup>t</sup>* � *<sup>k</sup> <sup>∂</sup>*2*<sup>u</sup> <sup>∂</sup>x*<sup>2</sup> ¼ *f x*ð Þ , *t* with *f x*ð Þ , *t* 6¼ 0 applies. 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 *u*ð Þ¼ 0, *t u L*ð Þ¼ , *t* 0. Given some initial profile of the temperature at time *t* ¼ 0, *u x*ð Þ¼ , 0 *h x*ð Þ, the analytical solution in 1D or 2D with a regular geometrical domain is given by separation of variables. Assuming *u x*ð Þ¼ , *t F x*ð Þ*G t*ð Þ, one obtains a system of two ODEs and from that, the series solution *u x*ð Þ¼ , *<sup>t</sup>* <sup>P</sup><sup>∞</sup> *<sup>n</sup>*¼<sup>1</sup>*un*ð Þ¼ *<sup>x</sup>*, *<sup>t</sup>* P<sup>∞</sup> *<sup>n</sup>*¼<sup>1</sup>*An* exp �*λ*<sup>2</sup> *nt* � � sin *<sup>n</sup><sup>π</sup> <sup>L</sup> <sup>x</sup>* � � with the coefficients *An* satisfying the initial condition *u x*ð Þ¼ , 0 <sup>P</sup><sup>∞</sup> *<sup>n</sup>*¼<sup>1</sup>*An* sin *<sup>n</sup><sup>π</sup> <sup>L</sup> <sup>x</sup>* � � <sup>¼</sup> *f x*ð Þ and the eigenvalues *<sup>λ</sup><sup>n</sup>* <sup>¼</sup> *nk<sup>π</sup> <sup>L</sup>* . The coefficients are evaluated via the sin series integration: *An* <sup>¼</sup> <sup>2</sup> *L* Ð *L* <sup>0</sup> *f x*ð Þsin *<sup>k</sup>π<sup>x</sup> L* � �*dx*. For nonintregrable 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].

**1.2 Numerical schemes**

*DOI: http://dx.doi.org/10.5772/intechopen.93928*

**2. Forward problem**

value problem in 1D:

operators:

**181**

the superposition of solutions of the form:

A quick calculation shows that *<sup>∂</sup>*

*∂*2 *u ∂x*2fd *∂u ∂t*fd

*xi*, *tj*

*xi*, *tj*

� � <sup>¼</sup> *u xi*þ1, *tj*

which results in the explicit finite difference scheme:

*ui*,*j*þ<sup>1</sup> ¼ *ui*,*<sup>j</sup>* þ

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 *u* at the next time step is given by a direct calculation from the previous one) or an implicit method (where the update for *u* is accomplished via

*Introduction to Numerical Approaches for Forward and Inverse Heat Transfer Problems*

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* ¼ 0 or in the inverse case, predict the temperature distribution at time *t* ¼ 0 from 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

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

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

�*π*2*kt* sin ð Þ *<sup>π</sup><sup>x</sup>*

� � � *u xi*, *tj*

*<sup>h</sup>*<sup>2</sup> *ui*þ1,*<sup>j</sup>* � <sup>2</sup>*ui*,*<sup>j</sup>* <sup>þ</sup> *ui*�1,*<sup>j</sup>*

*u x* ~ð Þ¼ , *t e*

*<sup>∂</sup><sup>t</sup>* � *<sup>k</sup> <sup>∂</sup>*<sup>2</sup> *∂x*<sup>2</sup>

� � <sup>¼</sup> *u xi*, *<sup>t</sup> <sup>j</sup>*þ<sup>1</sup>

*kl*

known explicit scheme, we use the forward in time and central in space difference

� � � <sup>2</sup>*u xi*, *tj*

*ut* ¼ *kuxx*, 0<*x*<1, (1) *u x*ð Þ¼ , 0 *f x*ð Þ, 0<*x*<1, (2) *u*ð Þ¼ 0, *t u*ð Þ¼ 1, *t* 0, *t*>0 (3)

h i*u x* <sup>~</sup>ð Þ¼ , *<sup>t</sup>* 0. In order to form the well-

*<sup>l</sup>* (4)

� � *<sup>h</sup>*<sup>2</sup> (5)

� � (6)

� �

� � <sup>þ</sup> *u xi*�1, *tj*

*<sup>∂</sup>x*<sup>2</sup> <sup>þ</sup> *<sup>∂</sup>*2*<sup>u</sup> ∂y*<sup>2</sup> � � <sup>¼</sup> *<sup>∂</sup><sup>u</sup>*

*<sup>∂</sup><sup>t</sup>* and

a linear solve). For 2D heat transfer, the model becomes *α*<sup>2</sup> *<sup>∂</sup>*2*<sup>u</sup>*

values. Regularization techniques need to be used for the solution [4].

**Figure 1.**

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

In this chapter we are concerned mainly with heat conduction for which the differential heat model applies. We take *Q* to 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: *<sup>Q</sup>*in � *<sup>Q</sup>*out <sup>¼</sup> *mc dT dt* , where *T* represents the temperature, *m* the mass, and *c* is the specific heat value of the material. It should be noted that in heat problems, *T* and *u* are 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: *<sup>q</sup>* ¼ �*<sup>k</sup> dT dx*, where *k* is the thermal conductivity of the material [2]. Using *Q*in ¼ *Q x*ð Þ and *Q*out ¼ *Q x*ð Þ þ *δx* along with heat flux *q* ¼ *Q=A*, we have *Aq x*ð Þ� *Aq x*ð Þ¼ <sup>þ</sup> *<sup>δ</sup><sup>x</sup> mc∂T=∂t*. Applying Fourier's law yields:

$$\overline{k}A\left[-\frac{\partial T}{\partial \mathbf{x}}\bigg|\_{\mathbf{x}} + \frac{\partial T}{\partial \mathbf{x}}\bigg|\_{\mathbf{x} + \delta \mathbf{x}}\right] = mc\partial T/\partial t$$

Next, using Taylor's approximation *<sup>∂</sup><sup>T</sup> ∂x <sup>x</sup>*þ*δx*≈*<sup>∂</sup><sup>T</sup> ∂x <sup>x</sup>* <sup>þ</sup> *<sup>∂</sup>*2*<sup>T</sup> ∂x*<sup>2</sup> *x δx* and expressing *m* ¼ *ρAδx* in terms of slab density, area, and width, we obtain:

$$
\overline{k}A\frac{\partial^2 T}{\partial \mathbf{x}^2} = \rho c \frac{\partial T}{\partial t} \quad \Rightarrow \quad k\frac{\partial^2 T}{\partial \mathbf{x}^2} = \frac{\partial T}{\partial 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=m*<sup>3</sup> (energy per unit volume). In this case, the in and out energy balance gives *Q x*ð Þþ *qV*\_ � *Q x*ð Þ¼ <sup>þ</sup> *<sup>δ</sup><sup>x</sup> mc∂T=∂t*, which with *<sup>q</sup>* <sup>¼</sup> *<sup>Q</sup>=A*,*<sup>V</sup>* <sup>¼</sup> *<sup>A</sup>δx*, Fourier's law and Taylor's expansion, gives *kA <sup>∂</sup>*2*<sup>T</sup> <sup>∂</sup>x*<sup>2</sup> *δx* þ *qA*\_ *δx* ¼ *mc <sup>∂</sup><sup>T</sup> <sup>∂</sup><sup>t</sup>* ) *<sup>k</sup> <sup>∂</sup>*2*<sup>T</sup> <sup>∂</sup>x*<sup>2</sup> <sup>þ</sup> *<sup>q</sup>*\_ <sup>¼</sup> *<sup>ρ</sup><sup>c</sup> <sup>∂</sup><sup>T</sup> <sup>∂</sup><sup>t</sup>* ) *<sup>k</sup> <sup>∂</sup>*2*<sup>T</sup> <sup>∂</sup>x*<sup>2</sup> <sup>þ</sup> *<sup>q</sup>*\_ *<sup>ρ</sup><sup>c</sup>* <sup>¼</sup> *<sup>∂</sup><sup>T</sup> <sup>∂</sup><sup>t</sup>* . In the case of steady state, *<sup>∂</sup>T=∂<sup>t</sup>* <sup>¼</sup> 0 and a second order differential equation model results.

**Figure 2.** *Heat conducting slabs without and with internal heat generation.* *Introduction to Numerical Approaches for Forward and Inverse Heat Transfer Problems DOI: http://dx.doi.org/10.5772/intechopen.93928*

### **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 *u* at the next time step is given by a direct calculation from the previous one) or an implicit method (where the update for *u* is accomplished via a linear solve). For 2D heat transfer, the model becomes *α*<sup>2</sup> *<sup>∂</sup>*2*<sup>u</sup> <sup>∂</sup>x*<sup>2</sup> <sup>þ</sup> *<sup>∂</sup>*2*<sup>u</sup> ∂y*<sup>2</sup> � � <sup>¼</sup> *<sup>∂</sup><sup>u</sup> <sup>∂</sup><sup>t</sup>* and 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* ¼ 0 or in the inverse case, predict the temperature distribution at time *t* ¼ 0 from 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**

In this chapter we are concerned mainly with heat conduction for which the differential heat model applies. We take *Q* to 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

*Perfectly insulated rod used a model of heat transfer, an exothermic slab, and sample modeled heat*

ture, *m* the mass, and *c* is the specific heat value of the material. It should be noted that in heat problems, *T* and *u* are ofter used interchangeably. We define *q* ¼ *Q=A* (heat over area) as the heat flux, which by Fourier's law is proportional to the

[2]. Using *Q*in ¼ *Q x*ð Þ and *Q*out ¼ *Q x*ð Þ þ *δx* along with heat flux *q* ¼ *Q=A*, we have

*∂x <sup>x</sup>*þ*δx*≈*<sup>∂</sup><sup>T</sup> ∂x <sup>x</sup>* <sup>þ</sup> *<sup>∂</sup>*2*<sup>T</sup> ∂x*<sup>2</sup> *x*

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=m*<sup>3</sup> (energy per unit volume). In this case, the in and out energy balance gives *Q x*ð Þþ *qV*\_ � *Q x*ð Þ¼ <sup>þ</sup> *<sup>δ</sup><sup>x</sup> mc∂T=∂t*, which with *<sup>q</sup>* <sup>¼</sup>

*<sup>∂</sup>x*<sup>2</sup> <sup>þ</sup> *<sup>q</sup>*\_

*<sup>ρ</sup><sup>c</sup>* <sup>¼</sup> *<sup>∂</sup><sup>T</sup>*

) *k*

*∂*2 *T <sup>∂</sup>x*<sup>2</sup> <sup>¼</sup> *<sup>∂</sup><sup>T</sup> ∂t*

*∂T ∂t*

*dt* , where *T* represents the tempera-

*δx* and expressing *m* ¼

*<sup>∂</sup>x*<sup>2</sup> *δx* þ *qA*\_ *δx* ¼

*<sup>∂</sup><sup>t</sup>* . In the case of steady state,

*dx*, where *k* is the thermal conductivity of the material

<sup>¼</sup> *mc∂T=∂<sup>t</sup>*

accordance with the model: *<sup>Q</sup>*in � *<sup>Q</sup>*out <sup>¼</sup> *mc dT*

*Inverse Heat Conduction and Heat Exchangers*

Next, using Taylor's approximation *<sup>∂</sup><sup>T</sup>*

*<sup>∂</sup>x*<sup>2</sup> <sup>þ</sup> *<sup>q</sup>*\_ <sup>¼</sup> *<sup>ρ</sup><sup>c</sup> <sup>∂</sup><sup>T</sup>*

*Heat conducting slabs without and with internal heat generation.*

*Aq x*ð Þ� *Aq x*ð Þ¼ <sup>þ</sup> *<sup>δ</sup><sup>x</sup> mc∂T=∂t*. Applying Fourier's law yields:

*kA* �*∂<sup>T</sup> ∂x x* þ *∂T ∂x x*þ*δx*

*ρAδx* in terms of slab density, area, and width, we obtain:

*<sup>Q</sup>=A*,*<sup>V</sup>* <sup>¼</sup> *<sup>A</sup>δx*, Fourier's law and Taylor's expansion, gives *kA <sup>∂</sup>*2*<sup>T</sup>*

*<sup>∂</sup>T=∂<sup>t</sup>* <sup>¼</sup> 0 and a second order differential equation model results.

*<sup>∂</sup><sup>t</sup>* ) *<sup>k</sup> <sup>∂</sup>*2*<sup>T</sup>*

*kA <sup>∂</sup>*<sup>2</sup> *T <sup>∂</sup>x*<sup>2</sup> <sup>¼</sup> *<sup>ρ</sup><sup>c</sup>*

temperature gradient: *<sup>q</sup>* ¼ �*<sup>k</sup> dT*

*mc <sup>∂</sup><sup>T</sup>*

**Figure 2.**

**180**

**Figure 1.**

*profile in 2D.*

*<sup>∂</sup><sup>t</sup>* ) *<sup>k</sup> <sup>∂</sup>*2*<sup>T</sup>*

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:

$$u\_t = k u\_{\infty}, \quad 0 < \infty < 1,\tag{1}$$

$$u(\mathbf{x}, \mathbf{0}) = f(\mathbf{x}), \quad \mathbf{0} < \mathbf{x} < \mathbf{1}, \tag{2}$$

$$u(\mathbf{0}, t) = u(\mathbf{1}, t) = \mathbf{0}, \quad t > \mathbf{0} \tag{3}$$

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:

$$
\tilde{u}(\mathbf{x},t) = e^{-\pi^2 kt} \sin\left(\pi \mathbf{x}\right),
$$

A quick calculation shows that *<sup>∂</sup> <sup>∂</sup><sup>t</sup>* � *<sup>k</sup> <sup>∂</sup>*<sup>2</sup> *∂x*<sup>2</sup> h i*u x* <sup>~</sup>ð Þ¼ , *<sup>t</sup>* 0. In order to form the wellknown explicit scheme, we use the forward in time and central in space difference operators:

$$\frac{\partial u}{\partial \mathbf{t}\_{\rm fd}} \left( \mathbf{x}\_{i}, t\_{j} \right) = \frac{u \left( \mathbf{x}\_{i}, t\_{j+1} \right) - u \left( \mathbf{x}\_{i}, t\_{j} \right)}{l} \tag{4}$$

$$\frac{\partial^2 u}{\partial \mathbf{x}^2 \epsilon \mathbf{j}} \left( \mathbf{x}\_i, t\_j \right) = \frac{u \left( \mathbf{x}\_{i+1}, t\_j \right) - 2u \left( \mathbf{x}\_i, t\_j \right) + u \left( \mathbf{x}\_{i-1}, t\_j \right)}{h^2} \tag{5}$$

which results in the explicit finite difference scheme:

$$u\_{i,j+1} = u\_{i,j} + \frac{kl}{h^2} (u\_{i+1,j} - 2u\_{i,j} + u\_{i-1,j}) \tag{6}$$

This can be re-written as *ui*,*j*þ<sup>1</sup> ¼ ð Þ 1 � 2*α ui*,*<sup>j</sup>* þ *α ui*þ1,*<sup>j</sup>* þ *ui*�1,*<sup>j</sup>* � � with *<sup>α</sup>* <sup>¼</sup> *kl <sup>h</sup>*2, with *l* the time step (Δ*t*) and *h* the spatial step size (Δ*x*). Von Neumann stability analysis yields the condition *l*≤ *<sup>h</sup>*<sup>2</sup> <sup>2</sup>*k*. For non-homogeneous cases the general model is *<sup>∂</sup><sup>u</sup> ∂t* ¼ *k <sup>∂</sup>*2*<sup>u</sup> <sup>∂</sup>x*<sup>2</sup> þ *f x*ð Þ , *t* and the finite difference analysis proceeds in a similar way, with an explicit scheme of the form: *ui*,*j*þ<sup>1</sup> <sup>¼</sup> *ui*,*<sup>j</sup>* <sup>þ</sup> *kl = h*2 *ui*þ1,*<sup>j</sup>* � 2*ui*,*<sup>j</sup>* þ *ui*�1,*<sup>j</sup>* � � <sup>þ</sup> *lf xi*, *tj* � �.

*d* 2

¼ *e*

¼ *e*

Setting *dt* � *kl*

obtain:

**183**

*u x* ~ð Þ , *t e*

<sup>¼</sup> <sup>2</sup>*kl*

1 2

nsfd scheme when *l* ≤ *<sup>h</sup>*<sup>2</sup>

**2.1 Numerical experiments**

�*π*2*kt* 1 <sup>2</sup>*<sup>i</sup> <sup>e</sup>*

> �*π*2*kt* 1 <sup>2</sup>*<sup>i</sup> <sup>e</sup> <sup>i</sup><sup>π</sup>xe <sup>i</sup>π<sup>h</sup>* � *<sup>e</sup>*

�*π*2*kt* 1 <sup>2</sup>*<sup>i</sup> <sup>e</sup> <sup>i</sup>π<sup>x</sup> e*

¼ *e*

*<sup>x</sup>*½ �¼ *u x* ~ð Þ , *t e*

*<sup>i</sup>π*ð Þ *<sup>x</sup>*þ*<sup>h</sup>* � *<sup>e</sup>*

*DOI: http://dx.doi.org/10.5772/intechopen.93928*

¼ *e*

¼ 2*e*

*h*2 *d* 2 *x*

�*π*2*kl* � <sup>1</sup> h i � <sup>2</sup>*kl*

*<sup>h</sup>*<sup>2</sup> ½ � cosð Þ� *<sup>π</sup><sup>h</sup>* <sup>1</sup>

*<sup>e</sup>*�*π*2*kl* � <sup>1</sup> cosð Þ� *πh* 1 " #

¼ *e*

�*iπ*ð Þ *<sup>x</sup>*þ*<sup>h</sup>* � � � 2 sin ð Þþ *<sup>π</sup><sup>x</sup>*

*Introduction to Numerical Approaches for Forward and Inverse Heat Transfer Problems*

<sup>2</sup>*<sup>i</sup>* 2 cosð Þ *<sup>π</sup><sup>h</sup> <sup>e</sup>*

�*i<sup>π</sup>xe*

*<sup>i</sup>π<sup>h</sup>* <sup>þ</sup> *<sup>e</sup>* �*iπ<sup>h</sup>* � � � *<sup>e</sup>*

�*π*2*kt* 1 2*i*

� �*u x* <sup>~</sup>ð Þ¼ , *<sup>t</sup>* 0, we obtain:

*kl <sup>h</sup>*<sup>2</sup> <sup>¼</sup> <sup>1</sup> 2

�*π*2*kt* 1

�*π*2*kt*½ � sin ð Þ� *<sup>π</sup>*ð Þ *<sup>x</sup>* <sup>þ</sup> *<sup>h</sup>* 2 sin ð Þþ *<sup>π</sup><sup>x</sup>* sin ð Þ *<sup>π</sup>*ð Þ *<sup>x</sup>* � *<sup>h</sup>* (10)

1 <sup>2</sup>*<sup>i</sup> <sup>e</sup>*

*<sup>i</sup>π<sup>h</sup>* � � (12)

�*iπ<sup>h</sup>* <sup>þ</sup> *<sup>e</sup> <sup>i</sup>π<sup>h</sup>* � � � <sup>4</sup>*i*sin ð Þ *<sup>π</sup><sup>x</sup>* � � (13)

�*π*2*kt* sin ð Þ *<sup>π</sup><sup>x</sup>* ½ �¼ cosð Þ� *<sup>π</sup><sup>h</sup>* <sup>1</sup> <sup>2</sup>*u x* <sup>~</sup>ð Þ , *<sup>t</sup>* ½ � cosð Þ� *<sup>π</sup><sup>h</sup>* <sup>1</sup> (16)

*<sup>h</sup>*<sup>2</sup> *u x* <sup>~</sup>ð Þ , *<sup>t</sup>* ½ �¼ cosð Þ� *<sup>π</sup><sup>h</sup>* <sup>1</sup> <sup>0</sup> ) *<sup>e</sup>*

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

*<sup>e</sup>*�*π*2*kl* � <sup>1</sup> cosð Þ� *πh* 1 " #

When *h*, *l* ≪ 1, the nsfd substitution coincides with the standard finite difference scheme. We use the small <sup>∣</sup>*q*<sup>∣</sup> Taylor series approximations *<sup>e</sup><sup>q</sup>*≈<sup>1</sup> <sup>þ</sup> *<sup>q</sup>* and sin ð Þ*<sup>q</sup>* <sup>≈</sup>*<sup>q</sup>* to

2

<sup>2</sup> just like in the original scheme. This also holds for large *h*,

<sup>1</sup> � *<sup>π</sup>*<sup>2</sup>*kl* � <sup>1</sup> � *<sup>π</sup>h*<sup>2</sup> 2

¼ *kl h*2

*<sup>h</sup>*<sup>2</sup> ≤1 is satisfied in the

" #

<sup>¼</sup> *<sup>e</sup>*�*π*2*kl* � <sup>1</sup> �2 sin <sup>2</sup> *<sup>π</sup><sup>h</sup>* 2 � � " #<sup>≈</sup> <sup>1</sup>

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, *u x*ð Þ¼ , 0 sin ð Þþ *πx x* for problem 1, with boundary conditions *u*ð Þ¼ 0, *t* 0 and *u*ð Þ¼ 1, *t* 1 and with the initial condition *u x*ð Þ¼ , 0 12 sin 9ð Þ� *πx* 7 sin 4ð Þ *πx* and zero homogeneous boundary conditions for problem 2. For the second example, errors vs. step size are also compared. In general, we

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

By this it follows that for small *h*, *l*, the stability constraint *kl*

as long as *h* is chosen so that ∣ cosð Þ� *πh* 1∣ is not very small.

observe more accurate results with the use of the nsfd scheme.

*<sup>i</sup><sup>π</sup>xe*

�*iπ<sup>x</sup>* � � � <sup>4</sup>*i*sin ð Þ *<sup>π</sup><sup>x</sup>* � � (14)

½ � 4*i* cosð Þ *πh* sin ð Þ� *πx* 4*i*sin ð Þ *πx* (15)

�*iπ*ð Þ *<sup>x</sup>*�*<sup>h</sup>* � � � � (11)

�*iπ<sup>h</sup>* � <sup>4</sup>*i*sin ð Þþ *<sup>π</sup><sup>x</sup> <sup>e</sup>*

�*iπ<sup>x</sup> e*

*<sup>i</sup>π<sup>x</sup>* � *<sup>e</sup>*

*<sup>i</sup>π*ð Þ *<sup>x</sup>*�*<sup>h</sup>* � *<sup>e</sup>*

�*iπ<sup>h</sup>* � *<sup>e</sup>*

�*i<sup>π</sup>xe*

�*π*2*kl* � <sup>1</sup> h i

If we use the backward time finite difference approximation for *<sup>∂</sup><sup>u</sup> ∂t* , we get: *ui*,*j*þ1�*ui*,*<sup>j</sup> <sup>k</sup>* <sup>¼</sup> *ui*þ1,*j*þ1�2*ui*,*j*þ1þ*ui*�1,*j*þ<sup>1</sup> *<sup>h</sup>*<sup>2</sup> , yielding the implicit scheme: �*αui*�1,*j*þ<sup>1</sup> <sup>þ</sup> ð Þ 1 þ 2*α ui*,*j*þ<sup>1</sup> � *αui*þ1,*j*þ<sup>1</sup> ¼ *ui*,*<sup>j</sup>* which can be solved for time *j* þ 1 from the information at time *j* using 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:

$$\frac{u\_{i,j+1} - u\_{i,j}}{l} = k \frac{\left(u\_{i+1,j+1} - 2u\_{i,j+1} + u\_{i-1,j+1}\right) + \left(u\_{i+1,j} - 2u\_{i,j} + u\_{i-1,j}\right)}{2h^2}$$

Again using *<sup>α</sup>* <sup>¼</sup> *kl <sup>h</sup>*2, we arrive at the scheme:

$$-a u\_{i+1,j+1} + \mathcal{D}(\mathbf{1} + a) u\_{i,j+1} - a u\_{i-1,j+1} = a u\_{i+1,j} + \mathcal{D}(\mathbf{1} - a) u\_{i,j} + a u\_{i-1,j}$$

Note that this takes the form of the linear system *Au* ¼ *b*, where *A* is a tridiagonal matrix of terms ð Þ �*α*,2 1ð Þ þ *α* , �*α* , as before, *u* containing the heat distribution at time *j* þ 1 and *b* containing 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:

$$\frac{T\_{i,j}^{n+1} - T\_{i,j}^n}{\Delta t} = k \left( \frac{T\_{i,j+1}^n - 2T\_{i,j}^n + T\_{i,j-1}^n}{\left(\Delta x\right)^2} + \frac{T\_{i+1,j}^n - 2T\_{i,j}^n + T\_{i-1,j}^n}{\left(\Delta y\right)^2} \right)$$

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:

$$
\overline{d}\_{\emptyset}f(\mathbf{x},t) = f(\mathbf{x}, t+l) - f(\mathbf{x}, t) \tag{7}
$$

$$
\overline{d}\_{\mathbf{x}}^{2}f(\mathbf{x},t) = f(\mathbf{x}+h,t) - 2f(\mathbf{x},t) + f(\mathbf{x}-h,t) \tag{8}
$$

and proceed to evaluate *dt* � *kl <sup>d</sup>* 2 *x h*2 � �*u x* <sup>~</sup>ð Þ , *<sup>t</sup>* for *u x* <sup>~</sup>ð Þ¼ , *<sup>t</sup> <sup>e</sup>*�*π*2*kt* sin ð Þ *<sup>π</sup><sup>x</sup>* . The stan-

dard 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 *kl <sup>h</sup>*2, which would force the equation to zero.

$$\begin{split} \overline{d}\_{l}[\tilde{u}(\mathbf{x},t)] &= \sin\left(\pi\mathbf{x}\right) \left[e^{-\pi^{2}k(t+l)} - e^{-\pi^{2}kl}\right] = \sin\left(\pi\mathbf{x}\right)e^{-\pi^{2}kl} \left(e^{-\pi^{2}kl} - \mathbf{1}\right) \\ &= \tilde{u}(\mathbf{x},t) \left(e^{-\pi^{2}kl} - \mathbf{1}\right) \end{split} \tag{9}$$

*Introduction to Numerical Approaches for Forward and Inverse Heat Transfer Problems DOI: http://dx.doi.org/10.5772/intechopen.93928*

$$\overline{d}\_{\mathbf{x}}^{2}[\tilde{\mu}(\mathbf{x},t)] = e^{-\pi^{2}kt} [\sin\left(\pi(\mathbf{x}+h)\right) - 2\sin\left(\pi\mathbf{x}\right) + \sin\left(\pi(\mathbf{x}-h)\right)] \tag{10}$$

$$=e^{-\pi^2kt}\left[\frac{1}{2i}\left(e^{i\pi(\mathbf{x}+h)}-e^{-i\pi(\mathbf{x}+h)}\right)-2\sin\left(\pi\mathbf{x}\right)+\frac{1}{2i}\left(e^{i\pi(\mathbf{x}-h)}-e^{-i\pi(\mathbf{x}-h)}\right)\right] \tag{11}$$

$$=e^{-\pi^2kt}\frac{1}{2i}\left[e^{i\pi\chi}e^{i\pi h}-e^{-i\pi\chi}e^{-i\pi h}-4i\sin\left(\pi\chi\right)+e^{i\pi\chi}e^{-i\pi h}-e^{-i\pi\chi}e^{i\pi h}\right] \tag{12}$$

$$=e^{-\pi^{2}kt}\frac{1}{2i}\left[e^{i\pi\chi}\left(e^{i\pi h}+e^{-i\pi h}\right)-e^{-i\pi\chi}\left(e^{-i\pi h}+e^{i\pi h}\right)-4i\sin\left(\pi\chi\right)\right] \tag{13}$$

$$=e^{-\pi^2kt}\frac{1}{2i}\left[2\cos\left(\pi h\right)\left(e^{i\pi\chi}-e^{-i\pi\chi}\right)-4i\sin\left(\pi\chi\right)\right] \tag{14}$$

$$=e^{-\pi^2kt}\frac{1}{2i}[4i\cos\left(\pi h\right)\sin\left(\pi\right)-4i\sin\left(\pi\infty\right)]\tag{15}$$

$$=2e^{-\pi^2kt}\sin\left(\pi x\right)[\cos\left(\pi h\right)-1]=2\tilde{u}(\pi,t)[\cos\left(\pi h\right)-1]\tag{16}$$

Setting *dt* � *kl h*2 *d* 2 *x* � �*u x* <sup>~</sup>ð Þ¼ , *<sup>t</sup>* 0, we obtain:

$$\begin{aligned} &\ddot{u}(\varkappa,t)\left[e^{-\varkappa^2kl}-\mathbf{1}\right]-\frac{2kl}{h^2}\ddot{u}(\varkappa,t)[\cos\left(\varkappa h\right)-\mathbf{1}]=\mathbf{0} \quad \Rightarrow \quad \left[e^{-\varkappa^2kl}-\mathbf{1}\right],\\ &\qquad =\frac{2kl}{h^2}\left[\cos\left(\varkappa h\right)-\mathbf{1}\right] \end{aligned}$$

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

$$\frac{kl}{h^2} = \frac{\mathbf{1}}{2} \left[ \frac{e^{-\pi^2 kl} - \mathbf{1}}{\cos\left(\pi h\right) - \mathbf{1}} \right]$$

When *h*, *l* ≪ 1, the nsfd substitution coincides with the standard finite difference scheme. We use the small <sup>∣</sup>*q*<sup>∣</sup> Taylor series approximations *<sup>e</sup><sup>q</sup>*≈<sup>1</sup> <sup>þ</sup> *<sup>q</sup>* and sin ð Þ*<sup>q</sup>* <sup>≈</sup>*<sup>q</sup>* to obtain:

$$\frac{1}{2} \left[ \frac{e^{-\pi^2 kl} - 1}{\cos \left( \pi h \right) - 1} \right] = \left[ \frac{e^{-\pi^2 kl} - 1}{-2 \sin^2 \left( \frac{\pi h}{2} \right)} \right] \approx \frac{1}{2} \left[ \frac{1 - \pi^2 kl - 1}{-\frac{\pi h^2}{2}} \right] = \frac{kl}{h^2}$$

By this it follows that for small *h*, *l*, the stability constraint *kl <sup>h</sup>*<sup>2</sup> ≤1 is satisfied in the nsfd scheme when *l* ≤ *<sup>h</sup>*<sup>2</sup> <sup>2</sup> just like in the original scheme. This also holds for large *h*, as long as *h* is 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, *u x*ð Þ¼ , 0 sin ð Þþ *πx x* for problem 1, with boundary conditions *u*ð Þ¼ 0, *t* 0 and *u*ð Þ¼ 1, *t* 1 and with the initial condition *u x*ð Þ¼ , 0 12 sin 9ð Þ� *πx* 7 sin 4ð Þ *πx* and 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

This can be re-written as *ui*,*j*þ<sup>1</sup> ¼ ð Þ 1 � 2*α ui*,*<sup>j</sup>* þ *α ui*þ1,*<sup>j</sup>* þ *ui*�1,*<sup>j</sup>*

If we use the backward time finite difference approximation for *<sup>∂</sup><sup>u</sup>*

tion at time *j* using a linear system with a tridiagonal matrix with elements

*<sup>l</sup>* <sup>¼</sup> *<sup>k</sup> ui*þ1,*j*þ<sup>1</sup> � <sup>2</sup>*ui*,*j*þ<sup>1</sup> <sup>þ</sup> *ui*�1,*j*þ<sup>1</sup>

*<sup>h</sup>*2, we arrive at the scheme:

yields the condition *l*≤ *<sup>h</sup>*<sup>2</sup>

*<sup>k</sup>* <sup>¼</sup> *ui*þ1,*j*þ1�2*ui*,*j*þ1þ*ui*�1,*j*þ<sup>1</sup>

*ui*,*j*þ<sup>1</sup> � *ui*,*<sup>j</sup>*

Again using *<sup>α</sup>* <sup>¼</sup> *kl*

*T<sup>n</sup>*þ<sup>1</sup> *<sup>i</sup>*,*<sup>j</sup>* � *<sup>T</sup><sup>n</sup> i*,*j* <sup>Δ</sup>*<sup>t</sup>* <sup>¼</sup> *<sup>k</sup>*

> *d* 2

and proceed to evaluate *dt* � *kl <sup>d</sup>*

*dt*½ �¼ *u x* ~ð Þ , *t* sin ð Þ *πx e*

¼ *u x* ~ð Þ , *t e*

implicit method.

equation to zero.

**182**

explicit scheme of the form: *ui*,*j*þ<sup>1</sup> <sup>¼</sup> *ui*,*<sup>j</sup>* <sup>þ</sup> *kl*

*Inverse Heat Conduction and Heat Exchangers*

*k <sup>∂</sup>*2*<sup>u</sup>*

*ui*,*j*þ1�*ui*,*<sup>j</sup>*

takes the form:

*l* the time step (Δ*t*) and *h* the spatial step size (Δ*x*). Von Neumann stability analysis

*<sup>∂</sup>x*<sup>2</sup> þ *f x*ð Þ , *t* and the finite difference analysis proceeds in a similar way, with an

*= h*2

*<sup>h</sup>*<sup>2</sup> , yielding the implicit scheme: �*αui*�1,*j*þ<sup>1</sup> <sup>þ</sup> ð Þ 1 þ 2*α ui*,*j*þ<sup>1</sup> � *αui*þ1,*j*þ<sup>1</sup> ¼ *ui*,*<sup>j</sup>* which can be solved for time *j* þ 1 from the informa-

ð Þ �*α*, 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

�*αui*þ1,*j*þ<sup>1</sup> þ 2 1ð Þ þ *α ui*,*j*þ<sup>1</sup> � *αui*�1,*j*þ<sup>1</sup> ¼ *αui*þ1,*<sup>j</sup>* þ 2 1ð Þ � *α ui*,*<sup>j</sup>* þ *αui*�1,*<sup>j</sup>*

*<sup>i</sup>*,*<sup>j</sup>* <sup>þ</sup> *<sup>T</sup><sup>n</sup> i*,*j*�1

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

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:

<sup>2</sup> þ

*Tn*

*dtf x*ð Þ¼ , *t f x*ð Þ� , *t* þ *l f x*ð Þ , *t* (7)

*xf x*ð Þ¼ , *t f x*ð Þ� þ *h*, *t* 2*f x*ð Þþ , *t f x*ð Þ � *h*, *t* (8)

¼ sin ð Þ *πx e*

!

*<sup>i</sup>*þ1,*<sup>j</sup>* � <sup>2</sup>*T<sup>n</sup>*

ð Þ Δ*y* 2

*u x* <sup>~</sup>ð Þ , *<sup>t</sup>* for *u x* <sup>~</sup>ð Þ¼ , *<sup>t</sup> <sup>e</sup>*�*π*2*kt* sin ð Þ *<sup>π</sup><sup>x</sup>* . The stan-

�*π*2*kt e*

*<sup>h</sup>*2, which would force the

(9)

�*π*2*kl* � <sup>1</sup> � �

*<sup>i</sup>*,*<sup>j</sup>* <sup>þ</sup> *<sup>T</sup><sup>n</sup> i*�1,*j*

Note that this takes the form of the linear system *Au* ¼ *b*, where *A* is a tridiagonal matrix of terms ð Þ �*α*,2 1ð Þ þ *α* , �*α* , as before, *u* containing the heat distribution at time *j* þ 1 and *b* containing 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:

ð Þ Δ*x*

2 *x h*2 � �

on the calculations below, we can make a replacement for *kl*

�*π*2*kl* � <sup>1</sup> � �

�*π*2*k t*ð Þ <sup>þ</sup>*<sup>l</sup>* � *<sup>e</sup>* �*π*2*kt* h i

dard 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

*<sup>i</sup>*,*j*þ<sup>1</sup> � <sup>2</sup>*T<sup>n</sup>*

*Tn*

<sup>2</sup>*k*. For non-homogeneous cases the general model is *<sup>∂</sup><sup>u</sup>*

� � <sup>þ</sup> *ui*þ1,*<sup>j</sup>* � <sup>2</sup>*ui*,*<sup>j</sup>* <sup>þ</sup> *ui*�1,*<sup>j</sup>*

2*h*<sup>2</sup>

*ui*þ1,*<sup>j</sup>* � 2*ui*,*<sup>j</sup>* þ *ui*�1,*<sup>j</sup>*

� � <sup>þ</sup> *lf xi*, *tj*

� �

� � with *<sup>α</sup>* <sup>¼</sup> *kl*

*<sup>h</sup>*2, with

*∂t* ¼

� �.

, we get:

*∂t*

In this case, separation of variables can be used and after a well-known sequence

� � exp �*<sup>k</sup> <sup>n</sup><sup>π</sup>*

*L* � �<sup>2</sup> *t*

� �. It follows that the sin coefficients of the solution at

� � which is a sin series; that is, *An* <sup>¼</sup>

*L* � �<sup>2</sup> *T* � �. The

*b u*<sup>0</sup> ð Þ,

ð Þ*t* ,

*unt* . The

*b* by the regularized

�ð Þ *<sup>j</sup>*þ<sup>1</sup> *<sup>U</sup>*<sup>0</sup> <sup>¼</sup> *<sup>u</sup> <sup>j</sup>*þ<sup>1</sup>

*u*2. This then implies that the initial

*<sup>i</sup>* � *<sup>u</sup> <sup>j</sup> i* <sup>Δ</sup>*<sup>t</sup>* � *Mu <sup>j</sup>*þ<sup>1</sup>

� �

*An* sin *<sup>n</sup>π<sup>x</sup> L*

*Introduction to Numerical Approaches for Forward and Inverse Heat Transfer Problems*

*L*

values of *u* at *t* ¼ *T* can 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

h i � � with the *:* <sup>∗</sup> representing component-wise opera-

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 *A* is a tridiagonal matrix of terms ð Þ �*α*,2 1ð Þ <sup>þ</sup> *<sup>α</sup>* , �*<sup>α</sup>* . This means that *<sup>u</sup>*<sup>1</sup> <sup>¼</sup> *<sup>A</sup>*�<sup>1</sup>

condition at time *<sup>t</sup>* <sup>¼</sup> 0, *<sup>u</sup>*<sup>0</sup> can be recovered from the application of the inverse matrix raised to a power *nt* corresponding to the time point. As another example, if

*<sup>i</sup>*�<sup>1</sup> ð Þ*<sup>t</sup>*

*<sup>h</sup>*<sup>2</sup> tridiagð Þ �1, 2, �<sup>1</sup> , for *<sup>j</sup>* <sup>¼</sup> 0, 1, … , *nt* � 1. Using *u x*ð Þ¼ , 0 *<sup>U</sup>*0, we

*<sup>u</sup> <sup>j</sup>* <sup>¼</sup> ð Þ *<sup>I</sup>* <sup>þ</sup> <sup>Δ</sup>*tM* �ð Þ *<sup>j</sup>*þ<sup>1</sup> *<sup>U</sup>*<sup>0</sup> ) *<sup>L</sup>*

� �, which for large enough *λ*>0, alleviates the

*<sup>i</sup>* ð Þþ*<sup>t</sup> <sup>u</sup> <sup>j</sup>*þ<sup>1</sup>

h i <sup>¼</sup> <sup>0</sup> <sup>¼</sup> *<sup>u</sup> <sup>j</sup>*þ<sup>1</sup>

This means that *<sup>U</sup>*<sup>0</sup> can in theory be obtained right from ð Þ *<sup>I</sup>* <sup>þ</sup> <sup>Δ</sup>*tM nt*

matrix needs not be formed explicitly, as the finite differencing scheme can be repeatedly applied to *<sup>u</sup>nt* <sup>¼</sup> *u T*ð Þ. 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 ð Þ *<sup>I</sup>* <sup>þ</sup> <sup>Δ</sup>*tM* �*nt* to the approximant for *<sup>b</sup>* <sup>¼</sup> *u T*ð Þ, 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

and utilize a regularization scheme.

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.

This can be accomplished by replacing the naive solution *L*�<sup>1</sup>

<sup>2</sup> <sup>þ</sup> *<sup>λ</sup>*∥*u*∥<sup>2</sup> 2

tion [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

*<sup>n</sup>*¼<sup>1</sup>*An* sin *<sup>n</sup>π<sup>x</sup>*

*<sup>t</sup>* <sup>¼</sup> *<sup>T</sup>* are given by the sin coefficients at *<sup>t</sup>* <sup>¼</sup> 0 multiplied by exp �*<sup>k</sup> <sup>n</sup><sup>π</sup>*

*t* ¼ *T*, *Fux* ½ � ð Þ , *T* , we can in principle go backwards, to obtain *u x*ð Þ¼ , 0

*b u*<sup>0</sup> ðÞ ) *b u*<sup>0</sup> ð Þ¼ *<sup>A</sup>*<sup>2</sup>

of steps, the following series solution can be derived [3]:

*n*¼1

*u x*ð Þ¼ , *<sup>t</sup>* <sup>X</sup><sup>∞</sup>

*L*

*L* � �<sup>2</sup> *T*

on quartile statistics of the nonzero absolute values.

*<sup>i</sup>*þ<sup>1</sup> ðÞ�*<sup>t</sup>* <sup>2</sup>*<sup>u</sup> <sup>j</sup>*þ<sup>1</sup>

It follows that *u x*ð Þ¼ , 0 <sup>P</sup><sup>∞</sup>

<sup>0</sup> *u x*ð Þ , 0 sin *<sup>n</sup>π<sup>x</sup>*

*DOI: http://dx.doi.org/10.5772/intechopen.93928*

2 *L* Ð *L*

*<sup>F</sup>*�<sup>1</sup> *Fux* ½ � ð Þ , *<sup>T</sup> :* <sup>∗</sup> exp *<sup>k</sup> <sup>n</sup><sup>π</sup>*

*b u*<sup>1</sup> ð Þ¼ *<sup>A</sup>*�<sup>2</sup>

*<sup>h</sup>*<sup>2</sup> *<sup>u</sup> <sup>j</sup>*þ<sup>1</sup>

*<sup>u</sup>*<sup>0</sup> <sup>¼</sup> *<sup>U</sup>*0; *<sup>u</sup> <sup>j</sup>*þ<sup>1</sup> <sup>¼</sup> ð Þ *<sup>I</sup>* <sup>þ</sup> <sup>Δ</sup>*tM* �<sup>1</sup>

�*nt*

we discretize *ut* ¼ *kuxx* as:

*<sup>∂</sup>u t*ð Þ *<sup>∂</sup><sup>t</sup>* � *<sup>k</sup>*

with *<sup>M</sup>* <sup>¼</sup> *<sup>k</sup>*

obtain the relations:

*LU*<sup>0</sup> ¼ *b* with *L* ¼ *L*

**185**

formulation min *<sup>u</sup>* <sup>∥</sup>*Lu* � *<sup>b</sup>*∥<sup>2</sup>

*Fux* ½ �¼ ð Þ<sup>0</sup>

*<sup>u</sup>*<sup>2</sup> <sup>¼</sup> *<sup>A</sup>*�<sup>1</sup>

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

#### **Figure 4.**

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

infinitesimal segment as in **Figure 1** yields the energy balance equation: *q*\_ þ *αdx* � �*<sup>q</sup>* <sup>þ</sup> *dq*\_ *dx dx* <sup>¼</sup> <sup>0</sup> ) *dq*\_ *dx* <sup>¼</sup> *<sup>α</sup>*. It follows that the heat model takes the form <sup>1</sup> *α ∂u <sup>∂</sup><sup>t</sup>* <sup>¼</sup> *<sup>∂</sup>*2*<sup>u</sup> <sup>∂</sup>x*<sup>2</sup> <sup>þ</sup> *<sup>q</sup>*\_ *k* . The steady state behavior does not depend on time so that *<sup>∂</sup><sup>u</sup> <sup>∂</sup><sup>t</sup>* ¼ 0 and an ordinary differential equation model *<sup>d</sup>*<sup>2</sup> *u dx*<sup>2</sup> <sup>þ</sup> *<sup>α</sup> <sup>k</sup>* ¼ 0 results. Upon a double integration, the general solution is *u x*ð Þ¼� *<sup>α</sup>* <sup>2</sup>*<sup>k</sup> <sup>x</sup>*<sup>2</sup> <sup>þ</sup> *Ax* <sup>þ</sup> *<sup>B</sup>*, which with the boundary conditions *<sup>T</sup>*ð Þ¼ <sup>0</sup> *T L*ð Þ¼ *Tw* becomes *u x*ð Þ¼� *<sup>α</sup>* <sup>2</sup>*<sup>k</sup> <sup>x</sup>*<sup>2</sup> <sup>þ</sup> *Lx* <sup>þ</sup> *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*∈½ � 0, *L* :

$$\frac{\partial u}{\partial t} = k \frac{\partial^2 u}{\partial \mathbf{x}^2} \tag{17}$$

$$
\mu(\mathbf{x}, \mathbf{0}) = f(\mathbf{x}), \quad \mu(\mathbf{0}, t) = \mathbf{0} = \mu(L, t) \tag{18}
$$

*Introduction to Numerical Approaches for Forward and Inverse Heat Transfer Problems DOI: http://dx.doi.org/10.5772/intechopen.93928*

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

$$u(\varkappa, t) = \sum\_{n=1}^{\infty} A\_n \sin\left(\frac{n\pi\varkappa}{L}\right) \exp\left(-k\left(\frac{n\pi}{L}\right)^2 t\right).$$

It follows that *u x*ð Þ¼ , 0 <sup>P</sup><sup>∞</sup> *<sup>n</sup>*¼<sup>1</sup>*An* sin *<sup>n</sup>π<sup>x</sup> L* � � which is a sin series; that is, *An* <sup>¼</sup> *Fux* ½ �¼ ð Þ<sup>0</sup> 2 *L* Ð *L* <sup>0</sup> *u x*ð Þ , 0 sin *<sup>n</sup>π<sup>x</sup> L* � �. It follows that the sin coefficients of the solution at *<sup>t</sup>* <sup>¼</sup> *<sup>T</sup>* are given by the sin coefficients at *<sup>t</sup>* <sup>¼</sup> 0 multiplied by exp �*<sup>k</sup> <sup>n</sup><sup>π</sup> L* � �<sup>2</sup> *T* � �. The values of *u* at *t* ¼ *T* can 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*, *Fux* ½ � ð Þ , *T* , we can in principle go backwards, to obtain *u x*ð Þ¼ , 0 *<sup>F</sup>*�<sup>1</sup> *Fux* ½ � ð Þ , *<sup>T</sup> :* <sup>∗</sup> exp *<sup>k</sup> <sup>n</sup><sup>π</sup> L* � �<sup>2</sup> *T* h i � � with the *:* <sup>∗</sup> 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 *A* is a tridiagonal matrix of terms ð Þ �*α*,2 1ð Þ <sup>þ</sup> *<sup>α</sup>* , �*<sup>α</sup>* . This means that *<sup>u</sup>*<sup>1</sup> <sup>¼</sup> *<sup>A</sup>*�<sup>1</sup> *b u*<sup>0</sup> ð Þ, *<sup>u</sup>*<sup>2</sup> <sup>¼</sup> *<sup>A</sup>*�<sup>1</sup> *b u*<sup>1</sup> ð Þ¼ *<sup>A</sup>*�<sup>2</sup> *b u*<sup>0</sup> ðÞ ) *b u*<sup>0</sup> ð Þ¼ *<sup>A</sup>*<sup>2</sup> *u*2. This then implies that the initial condition at time *<sup>t</sup>* <sup>¼</sup> 0, *<sup>u</sup>*<sup>0</sup> can be recovered from the application of the inverse matrix raised to a power *nt* corresponding to the time point. As another example, if we discretize *ut* ¼ *kuxx* as:

$$\frac{\partial u(t)}{\partial t} - \frac{k}{h^2} \left[ u\_{i+1}^{j+1}(t) - 2u\_i^{j+1}(t) + u\_{i-1}^{j+1}(t) \right] = 0 = \frac{u\_i^{j+1} - u\_i^j}{\Delta t} - M u^{j+1}(t),$$

with *<sup>M</sup>* <sup>¼</sup> *<sup>k</sup> <sup>h</sup>*<sup>2</sup> tridiagð Þ �1, 2, �<sup>1</sup> , for *<sup>j</sup>* <sup>¼</sup> 0, 1, … , *nt* � 1. Using *u x*ð Þ¼ , 0 *<sup>U</sup>*0, we obtain the relations:

$$\mathfrak{u}^0 = U\_0; \mathfrak{u}^{j+1} = \left(I + \Delta t \mathbf{M}\right)^{-1} \mathfrak{u}^j = \left(I + \Delta t \mathbf{M}\right)^{-\left(-j+1\right)} U\_0 \Rightarrow \overline{L}^{-\left(-j+1\right)} U\_0 = \mathfrak{u}^{j+1}$$

This means that *<sup>U</sup>*<sup>0</sup> can in theory be obtained right from ð Þ *<sup>I</sup>* <sup>þ</sup> <sup>Δ</sup>*tM nt unt* . The matrix needs not be formed explicitly, as the finite differencing scheme can be repeatedly applied to *<sup>u</sup>nt* <sup>¼</sup> *u T*ð Þ. 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 ð Þ *<sup>I</sup>* <sup>þ</sup> <sup>Δ</sup>*tM* �*nt* to the approximant for *<sup>b</sup>* <sup>¼</sup> *u T*ð Þ, 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 *LU*<sup>0</sup> ¼ *b* with *L* ¼ *L* �*nt* and utilize a regularization scheme.

This can be accomplished by replacing the naive solution *L*�<sup>1</sup> *b* by the regularized formulation min *<sup>u</sup>* <sup>∥</sup>*Lu* � *<sup>b</sup>*∥<sup>2</sup> <sup>2</sup> <sup>þ</sup> *<sup>λ</sup>*∥*u*∥<sup>2</sup> 2 � �, 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.

infinitesimal segment as in **Figure 1** yields the energy balance equation: *q*\_ þ *αdx* �

*Temperature distribution inside a thin slab with an exothermic material with homogeneous boundary*

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

steady state example with a finite difference scheme for the time derivative, yielding

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

> *∂*2 *u*

*u x*ð Þ¼ , 0 *f x*ð Þ, *u*ð Þ¼ 0, *t* 0 ¼ *u L*ð Þ , *t* (18)

The steady state behavior does not depend on time so that *<sup>∂</sup><sup>u</sup>*

*u dx*<sup>2</sup> <sup>þ</sup> *<sup>α</sup>*

1-D. Starting from the initial boundary value problem on *x*∈½ � 0, *L* :

*∂u <sup>∂</sup><sup>t</sup>* <sup>¼</sup> *<sup>k</sup>*

*dx* <sup>¼</sup> *<sup>α</sup>*. It follows that the heat model takes the form <sup>1</sup>

*α ∂u <sup>∂</sup><sup>t</sup>* <sup>¼</sup> *<sup>∂</sup>*2*<sup>u</sup> <sup>∂</sup>x*<sup>2</sup> <sup>þ</sup> *<sup>q</sup>*\_ *k* .

*<sup>∂</sup><sup>t</sup>* ¼ 0 and an ordinary

*<sup>k</sup>* ¼ 0 results. Upon a double integration, the gen-

*<sup>∂</sup>x*<sup>2</sup> (17)

<sup>2</sup>*<sup>k</sup> <sup>x</sup>*<sup>2</sup> <sup>þ</sup> *Ax* <sup>þ</sup> *<sup>B</sup>*, which with the boundary conditions *<sup>T</sup>*ð Þ¼ <sup>0</sup>

<sup>2</sup>*<sup>k</sup> <sup>x</sup>*<sup>2</sup> <sup>þ</sup> *Lx* <sup>þ</sup> *Tw*, a parabola. We can simulate this

�*<sup>q</sup>* <sup>þ</sup> *dq*\_ *dx dx* 

**184**

**Figure 4.**

**Figure 3.**

<sup>¼</sup> <sup>0</sup> ) *dq*\_

*conditions as simulated by an fd scheme.*

*Inverse Heat Conduction and Heat Exchangers*

differential equation model *<sup>d</sup>*<sup>2</sup>

*T L*ð Þ¼ *Tw* becomes *u x*ð Þ¼� *<sup>α</sup>*

the temperature profile shown in **Figure 4**.

eral solution is *u x*ð Þ¼� *<sup>α</sup>*

**3. Inverse problem**

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

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*∥<sup>∞</sup> and then iteratively lowered to a lower fraction (with the values either linearly or logarithmically spaced). At each *λ* value, the solution to *LTL* <sup>þ</sup> *<sup>λ</sup><sup>I</sup> <sup>u</sup><sup>λ</sup>* <sup>¼</sup> *Lt b* is 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<sup>λ</sup>* � *b*∥2*=N* reaches a plateau. It is possible to employ a *predictor-corrector* scheme to iteratively improve the solution. In this approach, the approximation to *U*<sup>0</sup> can be obtained as above, then the forward operator can be applied to propagate this approximated initial heat profile at *t* ¼ 0 to time *t* ¼ *T* and compare with observations [4]. We can compare the residual at time *T* and attempt some corrections to *U*<sup>0</sup> based on this residual such that the propagated solution at *t* ¼ *T* best matches observations. In the case of only a finite subset of observations, interpolation techniques can be utilized to complete the profile at *t* ¼ *T* across all values of *x* [7] prior to regularization to obtain an approximate smooth and continuous *u T*ð Þ heat profile.

½ � *N=*4, 3*N=*4 . For this, we take the sin transform of the solution at *t* ¼ *T* and

*Introduction to Numerical Approaches for Forward and Inverse Heat Transfer Problems*

*Initial condition reconstruction with different filtering levels and illustration of inversion for initial condition*

. In Octave, this can be accomplished with the code:

The term *R* may 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 *R* with either zero or some max value; e.g. taking the cutoff to be all values above 1*e*6 and replacing them either with 0 or 1*e*2. The cutoff can be estimated with standard statistical outlier techniques, by taking the quartiles *Q*1, *Q*2, *Q*<sup>3</sup> of the ∣*R*∣ values, computing *IQR* ¼ *Q*<sup>3</sup> � *Q*<sup>1</sup> and filtering all entries of *R* above *Q*<sup>3</sup> þ *αIQR* with *α*≥2 by setting them to the *Q*<sup>2</sup> value. The process is illustrated in **Figure 7**, for the same example as above with filtering for absolute value of entries above 1*e*1, 1*e*3, 1*e*5. It can be observed that filtering only the entries of *R* above 1*e*5 in 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 *u T*ð Þ both 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

This article discusses numerical approaches for the heat equation, both for the

homogeneous cases. For the forward heat propagation modeling problem, we have

forward and inverse problems with examples for homogeneous and non-

multiply by exp *k <sup>n</sup><sup>π</sup>*

R=**exp** (k.\*k <sup>0</sup>

**filter** (R);

**Figure 7.**

*t* ¼ 0 and *t* ¼ *T*.

**4. Conclusions**

**187**

*L* <sup>2</sup> *T*

*DOI: http://dx.doi.org/10.5772/intechopen.93928*

init\_cond\_approx=**real** (**ifft** (fuT.\*R));

\*T);

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

*using Tikhonov regularization with parameter tuned via the L-curve approach.*

#### **3.1 Numerical experiments**

As an example, consider the 1D cosine based initial heat profile with *k* ¼ 1 and 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 *u x*ð Þ¼ , 0 0*:*5 cos 2ð Þþ *π=*65*x* 1Þ, with support on

**Figure 6.** *Homogeneous 1D profile and change of temperature with time.*

*Introduction to Numerical Approaches for Forward and Inverse Heat Transfer Problems DOI: http://dx.doi.org/10.5772/intechopen.93928*

#### **Figure 7.**

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*∥<sup>∞</sup> and then iteratively lowered to a lower fraction (with the values either linearly or logarithmically spaced). At each *λ* value, the

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

As an example, consider the 1D cosine based initial heat profile with *k* ¼ 1 and 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 *u x*ð Þ¼ , 0 0*:*5 cos 2ð Þþ *π=*65*x* 1Þ, with support on

using the solution at the previously used *λ* as the initial guess. Typically, the progression to lower *λ* stops once the value of ∥*Lu<sup>λ</sup>* � *b*∥2*=N* reaches a plateau. It is possible to employ a *predictor-corrector* scheme to iteratively improve the solution. In this approach, the approximation to *U*<sup>0</sup> can be obtained as above, then the forward operator can be applied to propagate this approximated initial heat profile at *t* ¼ 0 to time *t* ¼ *T* and compare with observations [4]. We can compare the residual at time *T* and attempt some corrections to *U*<sup>0</sup> based on this residual such that the propagated solution at *t* ¼ *T* best matches observations. In the case of only a finite subset of observations, interpolation techniques can be utilized to complete the profile at *t* ¼ *T* across all values of *x* [7] prior to regularization to obtain an

approximate smooth and continuous *u T*ð Þ heat profile.

*Homogeneous 1D profile and change of temperature with time.*

*b* is obtained (e.g. with a Conjugate Gradient scheme)

solution to *LTL* <sup>þ</sup> *<sup>λ</sup><sup>I</sup> <sup>u</sup><sup>λ</sup>* <sup>¼</sup> *Lt*

*Inverse Heat Conduction and Heat Exchangers*

**Figure 5.**

**3.1 Numerical experiments**

**Figure 6.**

**186**

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

½ � *N=*4, 3*N=*4 . For this, we take the sin transform of the solution at *t* ¼ *T* and multiply by exp *k <sup>n</sup><sup>π</sup> L* <sup>2</sup> *T* . In Octave, this can be accomplished with the code:

fuT=**fft** (u (:, n)); *% temperature profile at time T* R=**exp** (k.\*k <sup>0</sup> \*T); **filter** (R); init\_cond\_approx=**real** (**ifft** (fuT.\*R));

The term *R* may 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 *R* with either zero or some max value; e.g. taking the cutoff to be all values above 1*e*6 and replacing them either with 0 or 1*e*2. The cutoff can be estimated with standard statistical outlier techniques, by taking the quartiles *Q*1, *Q*2, *Q*<sup>3</sup> of the ∣*R*∣ values, computing *IQR* ¼ *Q*<sup>3</sup> � *Q*<sup>1</sup> and filtering all entries of *R* above *Q*<sup>3</sup> þ *αIQR* with *α*≥2 by setting them to the *Q*<sup>2</sup> value. The process is illustrated in **Figure 7**, for the same example as above with filtering for absolute value of entries above 1*e*1, 1*e*3, 1*e*5. It can be observed that filtering only the entries of *R* above 1*e*5 in 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 *u T*ð Þ both 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* ¼ 0 and *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 nonhomogeneous 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.

**References**

2017.

Press, 2016.

Media, 2010.

(2002): 823–847.

(2010): 801–813.

**189**

[1] Özişik, M. Necati, Helcio RB Orlande, Marcelo J. Colaço, and Renato M. Cotta. Finite difference methods in heat transfer. CRC press,

*DOI: http://dx.doi.org/10.5772/intechopen.93928*

*Introduction to Numerical Approaches for Forward and Inverse Heat Transfer Problems*

[2] Glicksman, Leon R., and John H. Lienhard. Modeling and Approximation in Heat Transfer. Cambridge University

[3] Taler, Jan, and Piotr Duda. Solving direct and inverse heat conduction problems. Springer Science & Business

[4] Muniz, Wagner Barbosa, Haroldo F. de Campos Velho, and Fernando Manuel Ramos. "A comparison of some inverse methods for estimating the initial condition of the heat equation." Journal of Computational and Applied Mathematics 103, no. 1 (1999): 145–163.

[5] Mickens, Ronald E. "Nonstandard finite difference schemes for differential equations." Journal of Difference Equations and Applications 8, no. 9

[6] Hansen, Per Christian, and Dianne Prost O'Leary. "The use of the L-curve in the regularization of discrete ill-posed problems." SIAM journal on scientific computing 14, no. 6 (1993): 1487–1503.

[7] Castro, L. P., Q. Chen, and S. Saitoh. "Source inversion of heat conduction from a finite number of observation data." Applicable Analysis 89, no. 6
