Estimation of surface temperatures from interior measurements using Tikhonov regularization

(1)

Contents lists available atScienceDirect

Results in Applied Mathematics

journal homepage:www.elsevier.com/locate/results-in-applied-mathematics

Estimation of surface temperatures from interior

measurements using Tikhonov regularization

Jean Pierre Ngendahayo

a

_{, Japhet Niyobuhungiro}

a,∗

, Fredrik Berntsson

b a_{Department of Mathematics, School of Science, College of Science and Technology, University of Rwanda, P.O. Box}

3900 Kigali, Rwanda

b_{Department of Mathematics, Linköping University, SE-581 83 Linköping, Sweden}

a r t i c l e i n f o Article history:

Received 25 June 2020

Received in revised form 26 December 2020 Accepted 29 December 2020

Available online 7 January 2021

MSC: 49N45 93B40 35Q79 Keywords: Inverse problem Heat conduction Heat treatment Multiple measurements Tikhonov regularization a b s t r a c t

In many industrial applications, such as in the steel industry, it is important to monitor the surface temperature, and also heating or cooling rates, during an industrial process. However, often the surface may be inaccessible for direct measurements. Instead we are restricted to using interior measurements to estimate the surface temperature, by solving an inverse heat conduction problem. The problem is well-known and has been studied by many authors.

In this paper a novel algorithm is proposed, which allows the use of multiple measurements taken inside the material. The procedure is based on formulating the problem as an operator equation and stability is achieved by Tikhonov regularization. By an iterative procedure, we can deal with non-linear problems where the material properties depend on the temperature.

Numerical experiments show that the method works well. We also apply the proposed method to a real industrial problem with measured data taken during an steel heating experiment.

1. Introduction

The problem of estimating the surface temperature is of great importance in different industrial applications in engineering and science. We can mention applications such as in quenching studies, internal combustion engines, nuclear reactors, hypersonic flight vehicles, solid rocket nozzles, etc, [1,2]. Due to the thermal conditions it is often the case that the surface itself is inaccessible for direct measurements. Instead one is restricted to using the temperature histories measured at locations in the interior of the material to estimate the surface thermal conditions [3,4]. This situation occurs for example in relation to metal quenching applications, see [3,5].

The problem of estimating the surface temperature using interior measurements is often referred to as the inverse heat conduction problem (IHCP), see [4,6,7], and is an active area of research in inverse problems. The problem is ill-posed [8,9] in the sense that small errors in the recorded measurements may be amplified during the computations and can cause large errors in the solution [10]. Thus, regularization is needed to stabilize the computations. Various methods, for solving the IHCP have been studied in literature. For instance, transforming the problem into first-kind Volterra integral equation and then applying a standard regularization technique [11,12]. Also several methods based on filtering the data, e.g. by Meyer wavelets or Fourier transforms [6,9] or by using Mollification [7], have been considered. There are also methods

∗

Corresponding author.

E-mail address: jniyobuhungiro@ur.ac.rw(J. Niyobuhungiro). https://doi.org/10.1016/j.rinam.2020.100140

(2)

based on writing the exact solution in terms of a rapidly convergent series [13], see also [1,2,10]. Both linear and non-linear problems have been considered.

In this paper, the main idea is to design a linear operator mapping the surface temperature onto the temperature history recorded by the measurement device. A similar approach was used in [4] to improve the accuracy of a shielded thermocouple. In this paper the method is further improved by considering the use of multiple sets of measurements recorded inside the material. By discretizing the problem using the Crank–Nicolson implicit scheme we write the inverse problem as an overdetermined linear system of equations. Additional stability is also achieved by the use of Tikhonov regularization. While the method, strictly, requires the problem to be linear we also demonstrate that non-linear problems can be solved by adding a simple iterative procedure.

The paper is organized as follows: In Section 2we present the mathematical formulation of the considered IHCP. We also reformulate the problem as a linear operator equation and demonstrate its ill-posedness. We also show how to apply Tikhonov regularization to the problem. In Section3we apply the method to a numerical test problem with a known solution and show that it is effective and accurate. Further, in Section4we investigate a real world application with measured data, where the surface heat-flux on a steel slab is computed. Finally, in Section5we discuss our results as well as some potential improvements.

2. Problem specification and analysis 2.1. Formulation of the problem

Let us consider the one dimensional heat conduction equation: The temperature T (x

,

t) satisfies

(

κ

Tx)x

=

cp

ρ

Tt

,

0

<

x

<

L

,

(2.1)

where

κ

is the thermal conductivity [W

/

K m],

ρ

is the density [kg

/

m3_{], and c}

pis the specific heat capacity [J

/

kg K ]. In

this section we only consider the case when the equation is linear, i.e. the coefficients

κ

,

ρ

, and cpdepends on x but not

on the temperature T (x

,

t). The surface is located at x

=

0 and f (t)

=

T (0

,

t) is the sought after temperature history.

A measurement device, a thermocouple, is located at x

=

L and we assume that T (x

,

L)

=

g(t) is known. We also need

the initial temperature for the material. Several options are available, e.g. steady state, i.e. (

κ

Tx)x

|

t=0

=

0, or constant,

i.e. T (x

,

0)

=

T0. In this work we use the assumption that the initial temperature is constant. We do this mostly for

simplicity but also because it is appropriate for the application presented in Section4.

The model problem under consideration is thus the following: Find the temperature T (x

,

t) such that

⎧

⎪

⎨

⎪

⎩

(

κ

Tx)x

=

ρ

cpTt

,

t

≥

0

,

0

≤

x

≤

L

,

T (x

,

0)

=

T0

,

0

≤

x

≤

L

,

T (0

,

t)

=

f (t)

,

t

>

0

,

T (L

,

t)

=

g(t)

,

t

>

0

,

(2.2)

where f (t) is the sought after surface temperature and g(t) is known data.

If both f (t) and g(t), and also T0, are specified then the problem(2.2)is well-posed and the solution T (x

,

t) can be

computed. However, in order to find the unknown f (t), additional data needs to be available. In this work we consider an additional measurement device, located at x

=

L1, 0

<

L1

<

L, where the temperature history g1(t) is recorded. Formally,

the problem(2.2)can be used to define a mapping f (t)

↦→

g1(t)

:=

T (L1

,

t). However, unless T0

=

0 and g(t)

=

0 the

mapping would be non-linear, more specifically affine.

In order to obtain a linear mapping we first consider the problem: Find

v

(x

,

t) satisfying

⎧

⎪

⎨

⎪

⎩

(κv

x

)

x

=

ρ

cp

v

t

,

t

>

0

,

0

≤

x

≤

L

,

v

(x

,

0)

=

T0

,

0

≤

x

≤

L

,

v

(0

,

t)

=

T0

,

t

>

0

,

v

(L

,

t)

=

g(t)

,

t

>

0

.

(2.3)

Next, let the function u

=

T

−

v

satisfy the problem

⎧

⎪

⎨

⎪

⎩

(κ

ux

)

x

=

ρ

cput

,

t

>

0

,

0

≤

x

≤

L

,

u(x

,

0)

=

0

,

0

≤

x

≤

L

,

u(0

,

t)

=

f (t)

−

T0

,

t

>

0

,

u(L

,

t)

=

0

,

t

>

0

.

(2.4)

Then T (x

,

t)

=

u(x

,

t)

+

v

(x

,

t) is the solution of(2.2). The solution u(x

,

t) depends linearly on u(0

,

t)

=

f (t)

−

T0, and also u(L1

,

t)

=

g1(t)

−

v

(L1

,

t). This procedure is equivalent to setting g(t)

=

0 and T0

=

0 in(2.2). Thus we are motivated to

make the following definition.

Definition 2.1. Let T0

=

0 and g(t)

=

0, and also let 0

<

L1

<

L. For any function f (t)

∈

C0(

[

0

, ∞

)) we define an

operator K1

:

C0(

[

0

, ∞

))

−→

C0(

[

0

, ∞

)) by

g1(t)

=

(K1f )(t)

=

T (L1

,

t)

,

(2.5)

where T (x

,

t) is the solution to(2.2).

(3)

Note that if multiple measurements gi(t)

=

T (Li

,

t), 0

<

Li

<

L, for i

=

1

,

2

, . . . ,

m, are available we obtain an equal

number of operators Ki. Also we remark that the problem of finding T (x

,

t), from data T (L1

,

t)

=

g1(t) and T (L

,

t)

=

g(t)

is well-posed in the region L1

<

x

<

L, but ill-posed in the region 0

<

x

<

L1. The problem of finding the surface

temperature is thus reformulated as an ill-posed operator equation of the type (K1f )(t)

=

g1(t).

Two useful properties of the operator K1are given below. The proofs are straightforward and are omitted.

Lemma 2.2 (Linearity). The operator K1, fromDefinition 2.1, is linear. That is, for all f(1)

,

f(2)

∈

C0

([

0

, ∞

)

, and scalars

α

1,

α

2, the following property holds

(

K1(

α

1f(1)

+

α

2f(2))

)

(t)

=

α

₁

(

K1f(1)

)

(t)

+

α

₂

(

K1f(2)

)

(t)

.

Lemma 2.3 (Shift Invariance). The operator K1, fromDefinition 2.1, is time shift invariant. That is, for any t0

>

0, we have

K1f (t

−

t0)

=

g1(t

−

t0).

2.2. Discretization

In the definition of the operator K1, see(2.5), we have assumed that the functions involved, e.g. f (t), g(t) and g1(t),

belong to C(0)(

[

0

, ∞

)). In practice the observed measured temperature data will be available on a discrete set of grid points

{

tj

}

, j

=

0

,

1

, . . . ,

N. The data functions are thus replaced by vectors, e.g.,

f

=

(f (t1)

,

f (t2)

, . . . ,

f (tN))T

∈

RN

,

(2.6) where f (t0)

=

0 is the initial condition which is left out from the vector. Similarly we introduce the grid

{

xi

}

, i

=

1

,

2

, . . . ,

n, and introduce the discrete solution Tij

=

T (xi

,

tj). In our experiments we assume that the grid is equidistant

in both x- and t-coordinates, but not necessarily with the same step size.

The boundary value problem (2.2)can be solved using a standard second order accurate finite difference scheme. For our experiments we use the Crank–Nicolson method [14]. Recall that we are interested in solving problems with temperature dependent material properties, e.g.

κ := κ

(T ). Since the solution depends on both x and t this is equivalent to the material properties also being a function of both x and t, e.g.

κ

(x

,

t)

:=

κ

(T (x

,

t)). In our finite difference code we

thus assume that the coefficients are given as matrices

{

κ

ij

}

,

{

ρ

ij

}

, and

{

(cp)ij

}

, for 1

≤

i

≤

n and 1

≤

j

≤

N. If the coefficients

are allowed to depend on both x and t the operator(2.5)is still linear but the time shift property, i.e.Lemma 2.3, does not hold.

The solution to the model problem (2.2) can be computed using Crank–Nicolson scheme even in the case of temperature dependent material properties. The non-linearity is simply dealt with by fixed point iteration where an initial guess T(0) is used to compute the initial approximation for the material properties, e.g.

κ

(0). The boundary value problem is then solved using the finite difference method giving us the next approximate solution T(1)_{. The process is}

rather efficient and converges rapidly [4].

In the discrete setting both f and g1are vectors in RN, cf.(2.6). Since the operator K1is linear, and also the numerical

solution, computed by Crank–Nicolson for fixed coefficient matrices

{

κ

_ij

}

,

{

ρ

_ij

}

, and

{

(cp)ij

}

, depend linearly on the data

vector f , we actually have a linear system of equations,

K1f

=

g1

,

K1

∈

RN×N

.

(2.7)

In this work, for simplicity, we will use the same symbol to denote the operator and the matrix, and also the vectors and the respective functions.

In the case where the material properties depend on x but are independent of temperature we have the following structure for the matrix K1.

Lemma 2.4. The matrix K1in(2.7)is a lower triangular Toeplitz matrix of the form

K1

=

⎡

⎢

⎣

k1 0 0 0

· · ·

0 k2 k1 0 0

· · ·

0 k3 k2 k1 0

· · ·

0

...

... ... ... ...

kN−1 kN−2 kN−3

... ...

0 kN kN−1 kN−2

· · ·

k1

⎤

⎥

⎦

.

(2.8)

Proof. Let

{

ei

}

be the standard basis of the space RN. Using the linearity of K1we have

K1

(

_N

∑

i=1 f (ti)ei

)

=

N

∑

i=1 f (ti)K1ei

=

g1

.

(2.9) 3

(4)

Thus, in order to evaluate the matrix–vector product K1f we need the vectors K1ei, i

=

1

,

2

, . . . ,

N. Here K1e1

=

K1(

:

,

1) is

the first column of the matrix K1. The remaining columns are obtained using the shift invariance, i.e.Lemma 2.3. In fact,

the entries kjin the matrix K1in(2.8)are given by kj

=

T

(

L1

,

tj

)

, j

=

1

, . . . ,

N, where T is the numerical solution to the

boundary value problem(2.2)obtained via the Crank Nicholson scheme, with data f

=

e1, g

=

0 and T0

=

0. Thus the

matrix K1has the advertised structure. □

The important point in the above lemma is that in order to compute the entire matrix K1 it is sufficient to solve the

boundary value problem(2.2)once, with data f

=

e1, g

=

0 and T0

=

0.

2.3. Multiple measurement points

In Section 2.1we investigated the use of two thermocouples, located at x

=

L and x

=

L1

<

L, providing the data T (L1

,

t)

:=

g1(t) and T (L

,

t)

=

0. Thus, in the discrete setting, the unknown surface temperature can in principle be

obtained by solving the linear system K1f

=

g1.

In this section we make the observation that additional measurements, i.e. vectors gi

:=

T (Li

,

t), Li

<

L, for i

=

2

,

3

, . . . ,

m, might be available. If multiple thermocouple locations are available we obtain an equal number of matrices Ki, such that gi

=

Kif . The linear system K1f

=

g1is thus replaced by an overdetermined system of equations,

K(m)f

=

g(m)

,

where K(m)

=

⎛

⎜

⎝

K1 K2

...

Km

⎞

⎟

⎠

,

and g(m)

=

⎛

⎜

⎝

g1 g2

...

gm

⎞

⎟

⎠

,

(2.10)

and solve the corresponding least squares problem

min f∈RN m

∑

i=1

∥

Kif

−

gi

∥

22

.

(2.11)

It is reasonable to assume that by using multiple measurements in estimation of surface temperature we improve the accuracy of the solution. In particular in the presence of random measurement noise. Note also that all the matrices

Ki may be computed by solving the problem(2.2), with data f

=

e1, g

=

0 and T0

=

0, once. After all, Ki

(:,

1

) =

(

T (Li

,

t1)

, . . . ,

T (Li

,

tN)

)

T, for i

=

1

, . . . ,

m, where T is the numerical solution to the boundary value problem(2.2)obtained

via the Crank Nicholson scheme with data f

=

e1, g

=

0 and T0

=

0.

2.4. Ill-posedness and regularization

The problem of estimating the surface temperature of a body, given interior measurements, is ill-posed in the sense that small errors in the measurements can cause large errors in the computed solution [6]. As a consequence the discrete problem, i.e. the linear systems(2.7)or(2.10), might be severely ill-conditioned.

In practice any measurement will contain errors, i.e. we actually have data

g₁δ

=

g1

+

η,

(2.12)

where

η

is random noise, representing the error in the data, with known noise level

∥

g₁δ

−

g1

∥

2

≤

δ

. By making use of

the singular value decomposition of the matrix K1, we can write the numerical solution as

fδ

=

N

∑

i=1

v

T ig1δ

σ

i ui

=

N

∑

i=1

v

T ig1

σ

i ui

+

N

∑

i=1

v

T i

η

σ

i ui

,

(2.13)

where

{

ui

}

and

{

v

i

}

are the singular vectors and

{

σ

i

}

are the singular values. The second term in(2.13)represents the

influence of random noise on the solution. Here small singular values are problematic and since the original problem is ill–posed we expect the discrete problem to have very small singular values

σ

i. In fact, singular values in such problems

decay rapidly to 0. This is typical for discrete ill-posed problems [15].

Example 2.5. In order to investigate the singular values for the discrete matrix K(m)_{, see}_(2.10)_{, we select a uniform grid}

of size N

=

103_{on the time interval 0}

_<

_t

_<

_{5, and compute the matrices K}(m)_{by solving the boundary value problem}

(2.2)using our Crank–Nicolson code. For the space discretization a uniform grid of size n

=

200 was used and the second boundary was located at L

=

1. The material parameters were set to

ρ =

1, cp

=

1, and

κ

(x)

=

1

+

0

.

3 cos(x2)

,

0

<

x

<

1

,

(2.14) so that the equation is linear. We consider two cases, each with three additional thermocouples, namely

(L1

,

L2

,

L3)

=

(0

.

70

,

0

.

75

,

0

.

80) and (L1

,

L2

,

L3)

=

(0

.

70

,

0

.

90

,

0

.

95)

.

(2.15)

(5)

Fig. 1. We display the first 200 singular values for the matrices K(m)_{as detailed in}_{Example 2.5. In all cases we show the case m}₌_{1 (blue curves),}

m=2 (red curves) and m=3 (black curves). We also display the case when the thermocouple locations are Li=0.70,0.75, and 0.80 (left graphs) and the case when Li=0.70,0.90, and 0.95 (right graphs). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We compute the singular values of the matrices K(m)_{, for m}

₌

₁

_,

₂

_,

_{3, and for both sets of thermocouple locations. The}

results are shown inFig. 1. The same matrix K1

=

K(1) is used. However since the locations L2and L3differ, so will the

matrices K(2) and K(3). The results show that in all cases the rate of decay for the singular values is similar. However, the matrix K(3) _{has singular values of a larger magnitude than the matrix K}(1)_{. Thus the ill-conditioning of the problem is}

improved by using more measurement points. We also see from the results that it is better to add measurement points as close to the unknown surface as possible.

In order to stabilize the computations we use Tikhonov regularization [16], and add an extra term to the least squares problem(2.11). Thus the regularized solution is given by

f_λδ

=

arg min f∈RN

(

_m

∑

i=1

∥

Kif

−

giδ

∥

2 2

+

λ

2

_∥

_f

_∥

2 2

)

,

(2.16)

where the vectors g_iδ, i

=

1

,

2

, . . . ,

m represent the noisy measurements at location x

=

Li, and

λ >

0 is the regularization

parameter. By using Tikhonov regularization the singular values

{

σ

i

}

, of the matrix K(m)are effectively replaced by the

sequence

{

σ

_i(

σ

2

i

+

λ

2)

−1

_}

_{which is bounded. Thus the effects of random noise on the computations is limited. The}

appropriate choice of the regularization parameter

λ

represents the trade off between solving the problem accurately, i.e. keeping the residual small, and limiting the magnification of random noise. There are several good methods for choosing the regularization parameter and in this work we will use the L-curve [15].

3. Numerical experiments

In this section, we present some numerical experiments in order to illustrate the usefulness of using multiple measurement points in estimation of surface temperature. We use the Crank–Nicolson scheme to create numerical test problems with a known solution, and also to create the matrices Ki, i

=

1

,

2

, . . . ,

m, representing the linear operators

involved.

For all tests the measurements are assumed to take place on the interval 0

<

t

<

10, and the boundary conditions

T (x

,

0)

=

0 and T (L

,

t)

=

0, where L

=

2, are used, in the problem(2.2). Thus we avoid the issue of solving the auxiliary problem(2.3)in order to set g(t) and T0to zero. Also all the tests presented in this section use a grid size of N

=

1000

for the time variable, and a grid size of n

=

200 in space. Also, in this section we only investigate the linear case and set

ρ =

1, cp

=

1, and

κ

(x)

=

1

+

0

.

3 cos(x2₎

_,

₀

_<

_x

_<

₂

_.

_(3.17)

Example 3.1. In our first test we use a rather simple function f (t) as the exact solution and solve the problem(2.2)in

order to simulate measurement data. The exact solution f (t) and the simulated measurements gi(t), at locations x

=

Li,

for Li

=

1

.

0, 1.05 and 1.10, are displayed inFig. 2. We also show the noisy data vectors after adding normally distributed

noise of variance 2

·

10−3_.

First we test the Tikhonov method for the case of a single measurement point. In this case we use the data vector g₁δ, recorded at L1

=

1

.

00, and try to reconstruct the surface temperature f (t). In order to find the appropriate regularization

(6)

Fig. 2. We display the exact solution f (t) (black, dashed) and the simulated measurements g1(t) (blue, solid), g2(t) (red, solid) and g3(t) (black, solid). We display both the exact, noise-free, data (left graphs) and the data with normally distributed noise, of variance 2·10−3_{, added (right graphs). This} is a rather significant amount of noise for this test. The measurements are recorded at Li=1.00,1.05, and 1.10. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

parameter

λ

we compute the reconstruction error

∥

f_λδ

−

f

∥

₂, as a function of

λ

. The results show that

λ =

2

·

10−2_{gives a}

suitable level of regularization. We also display the L-curve and mark its corner. The regularization parameter, suggested by the L-curve, is

λ =

2

.

5

·

10−2. This is close to the optimal value. The results are shown inFig. 3.

Second, we use both two and three measurement locations. In the case of two locations, e.g. L1

=

1

.

00 and L2

=

1

.

05

the smallest error is achieved for

λ =

4

·

10−2_{, and if all three locations are used the optimal value is}

_{λ =}

₅

_·

₁₀−2_{. The results}

are displayed inFig. 4. We see that more measurements used leads to a smaller error. The improvement after adding the second measurement locations is much larger than the improvement form adding the third set of measurements.

Example 3.2. In our second experiment we use a discontinuous function as our exact solution f (t). We also select 9 equally

spaced measurement points Li in the interval 1

.

00

<

Li

<

1

.

08, and compute the simulated data giδ(t), i

=

1

,

2

, . . . ,

9.

The added, normally distributed noise had a variance 5

·

10−3_{. This is slightly larger than the previous experiment. The}

appropriate level of regularization is selected by computing the error

∥

f_λ(m)

−

f

∥

2, for a large range of parameter values

λ

, and for all possible choices for the number of measurement locations m. The optimal

λ

, and also the smallest error achieved, as a function of m is displayed inFig. 5. We see that for a larger m we also need a larger

λ

. This is natural since a single penalty term

λ

2

∥

f

∥

2₂has to balance several residual terms

∥

Kif

−

giδ

∥

2

2. We also see that by increasing m

we reduce the error by a fair amount. Adding measurements close to L1

=

1

.

00 gives a larger reduction than adding

measurements taken further away. We also display several solutions for different values of m. The results show that adding more measurements does indeed help improve the computed solution.

4. Heating rates during heat treatment of steel

Heat treatment of steel is an important industrial process that needs careful control in order for the final product to have the desired properties. The experiment presented in this Section was carried out using a laboratory scale furnace and the steel slab was of dimension 100

×

80

×

22 [mm]. Thermocouples were placed in the middle of the slab, i.e. at (50

,

40) mm, and at distance L1

=

5 mm, L2

=

11 mm and L3

=

17 mm from the surface of the slab. The thermocouples

are referred to as TC1, TC2and TC3respectively. During the experiment the slab, initially at room temperature, was heated

by a gas burner. The temperature at the measurement locations were recorded using a sampling rate of 1

/

10 Hz and the data vectors are of size N

=

656. Thus the experiment lasted a little less than two hours. Our aim is to reconstruct both the temperature and the heat flux at the surface.

The furnace and also the recorded data are shown inFig. 6. The steel is of a known composition, with known material properties in the temperature range relevant for the experiment. Both the thermal conductivity

κ

and also the product

ρ

cpare shown inFig. 7. Since the problem is non-linear the Tikhonov regularization, see(2.16), needs to be combined

with an iterative scheme where a linear problem is solved in each step. We have three measurement locations L1, L2, and L3. The location L3takes the role of L in the definition of the linear operator, see2.1, and thus the Tikhonov procedure,

see(2.16), will use m

=

2 linear operators K1and K2.

In this application the exact solution f (t)

=

T (0

,

t) should be smooth. However, the temperature range is quite large,

i.e. 20–800o_{C and thus the norm}

_∥

_f

_∥

2

2 may also be quite large. Since the Tikhonov regularization, see(2.16), uses the

penalty term

λ

2

_∥

_f

_∥

2

2, we tend to obtain solutions fλ(t) with too small temperatures. Thus the heat flux close to the surface

(7)

Fig. 3. We display the error∥f_λδ−f∥2as a function ofλ(top, left) for the case of using a single measurement point located at L1=1.00. We also show the L-curve (top, right) for the same case. The corner of the L-curve is located atλTik≈3.5·10−2. We also display the numerical solutions f_λδ(t) for the caseλ = λTik(bottom, left), which corresponds to the minimum error, and for the caseλ =3·10−3(bottom, right). In both cases the exact solution is shown (black) and the numerical solution (red). It is clear thatλ = λTikchosen by L-curve corresponds to the best stable solution (bottom, left). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. We display the exact solution f (t) (black) and the numerical solution f_λδ(t) (red) computed using three measurement points. For this case the regularization parameterλ =5·10−2_{was used. We also show the point wise errors}_{|f (t)}₋_fδ

λ(t)|for the case of m=1 measurement location andλ =2·10−2_{(blue), m}₌_{2 locations and}_{λ =}₄_·₁₀−2_{(red), and m}₌_{3 with}_{λ =}₅_·₁₀−2 _{(black). Clearly, we see that more measurements} used leads to a smaller error. The measurements were recorded at Li=1.00,1.05, and 1.10. Moreover, the computed solution (left) is clearly more stable than the computed solution inFig. 3(bottom, left). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

(8)

Fig. 5. The optimal regularization parameterλ(top, left) and the smallest achieved error∥f_λδ−f∥2 (top, right), as a function of the number of measurements used. We also display the solutions f_λδ(t) (red) and the exact solution f (t) (black), for the optimal value of the regularization parameter

λ, and for the case m=1 (middle, left), m=3 (middle, right), m=5 (bottom, left) and m=9 (bottom, right). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

is also underestimated. An alternative is to use a penalty term of the form

λ

2

_∥

_f′′

(t)

∥

2_{. Thus we introduce a matrix D} 2, which

represents a central difference approximation of the second derivative on the grid, and instead solve the minimization problem f_λδ

=

arg min f∈RN

(

_m

∑

i=1

∥

Kif

−

giδ

∥

2 2

+

λ

2

_∥

_D 2f

∥

22

)

.

(4.18) 8

(9)

Fig. 6. The furnace used in the experiment (left) and the recorded data (right). The data recorded by the thermocouple TC1(blue), TC2 (red) and

TC3(black) are shown. Note that the curves are smooth and relatively little random noise is present. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7. We also show the material propertiesκandρcp, as a function of the temperature, for the specific steel quality used.

This will force the solution f_λδ to be smooth but not limit its magnitude. This is also very effective in preventing random noise from entering into the calculations. We refer the reader to [17], for higher-order Tikhonov regularization.

For the computations we used a regularization parameter

λ =

7

·

10−3_{, which is sufficient to produce a very smooth}

surface temperature f_λ(t) but leaves a small random noise visible in the computed thermal gradients Tx(0

,

t). The

well-posed heat conduction problems were solved using Crank–Nicolson and a grid of size n

=

200 in the x-direction. Our procedure is detailed below and the first step is used to illustrate the numerical scheme.

The starting guess is obtained by observing that the temperature range in the experiment is roughly 20oC to 800oC. Thus we simply use T(0)(x

,

t)

=

400oC, as our initial guess and calculate the initial approximation of the coefficients

κ

(0)_(x

i

,

tj) and (

ρ

cp)(0)(xi

,

tj).

A general step is as follows: Given the coefficients

κ

(k)_(x

i

,

tj) and (

ρ

cp)(k)(xi

,

tj), we first estimate the initial temperature,

T0of the slab. We know that the slab is initially at a constant temperature and the thermocouples only record random

noise. In our computations we set T0to be the mean value of the first two measurements. Then we solve the auxiliary

problem(2.3), with data

v

(k)_(x

_,

₀₎

₌

_T

0,

v

(k)(0

,

t)

=

T0, and

v

(k)(L3

,

t)

=

TC3(t). Note that, since the coefficients are now

functions of (x

,

t), and not of the temperature T , the equation is linear. The solution

v

(0)_(x

_,

_{t) obtained during the first}

step of the procedure is displayed inFig. 8. The second step consists of computing the data for the auxiliary problem(2.4). Thus we compute new data functions g₁(k)(t)

=

TC1(t)

−

v

(k)(L1

,

t) and g2(k)(t)

=

TC2(t)

−

v

(k)(L2

,

t). We also calculate the

two matrices K₁(k)and K₂(k)at this point. Since the coefficients are of the form

κ

(x

,

t), the time shift property does not hold

(10)

Fig. 8. Illustration of the iterative procedure by displayingv(0)_(x_,_{t) (top-left) and the data g}(0)

1 (t) (blue, top-right) and g (0)

2 (t) (red, top-right). Also

u(0)_(L3_,_t) ₌_TC_3(t)₋_v(0)_(L3_,_{t) is shown (black,top-right). We also show the solution u}(0)_(x_,_{t) (bottom,left) and the next approximation of the} temperature T(1)_(x_,_{t) (bottom,right). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of} this article.)

and we need to compute each matrix–vector product K_ℓej, for

ℓ =

1

,

2, j

=

1

,

2

, . . . ,

N, separately. This can still be done

using a single run of the Crank–Nicolson scheme, if the code is written to allow for using multiple right hand sides. The updated data functions g₁(0)(t) and g₂(0)(t) are shown inFig. 8.

The third step is selecting a regularization parameter

λ

and finding the minimizer f_λ(k)(t), of the Tikhonov functional

(4.18). We then solve the heat conduction problem(2.4), with data u(0

,

t)

=

f_λ(k)(t), u(x

,

0)

=

0, and u(L3

,

t)

=

0, and

compute the next approximation of the temperature T(k+1)

=

u(k)

+

v

(k). We also compute new coefficients

κ

(k+1)and (

ρ

cp)(k+1)at this point. The computed solution u(0)(x

,

t) and the next temperature approximation T(1)(x

,

t) are also shown

inFig. 8.

Finally, we check for convergence by computing the Frobenius norm

∥

T(k+1)

−

T(k)

∥

F. We can then start the procedure

again for the next iteration.

In Fig. 9we show both the final solution f_λ(10)(t)

:=

T(10)₍₀

_,

_{t), obtained using}

_{λ =}

₁₀−2_{, and also the convergence}

history

∥

T(k+1)

₋

_T(k)

_∥

F during the iterations. We remark that since, both the solution and the data are smooth, does not

contain any sharp features, and also has relatively small amounts of noise, the solution f_λ(k) is not at all sensitive with respect to the choice of

λ

. Also the rate of convergence is rapid and just two or three iterations produce identical results. A Frobenius norm

∥

T(3)

₋

_T(2)

_∥

F

≈

25

.

2 means an average difference of about 1

.

9

·

10−4at the grid points. This is too

small to be noticed in graphs.

We also compute the heat flux

κ

Tx(x

,

t), for the surface x

=

0 and for x

=

L1, using a one sided finite difference

quotient. We see that, while the difference is not very big, the heat-flux at the surface has a sharper peak compared to the curve 5 mm below the surface. The results are displayed inFig. 9.

5. Concluding remarks

In this paper we have considered the well-known problem of estimating the surface temperature using measurements recorded inside a material. The problem is frequently treated as a Cauchy problem [9], or as an operator equation [12], and several good numerical schemes exist. The novel feature of our proposed method is that it allows the use of multiple

(11)

Fig. 9. The computed solution f_λ(10)(t) (red, top-left), usingλ =10−2_{, and also the measurements TC1(t) (black, top-left) are displayed. Also the} convergence history (top-right) during the iterations is presented. Further, we display the computed heat-flux at the surface , i.e.κTx(0,t) (bottom-left), and at the location of the first thermocouple, i.e.κTx(L1,t) (bottom-right). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

measurements, recorded at different locations, inside the material. This allows us to estimate the surface temperature with a much lower error. In its basic form the method requires the heat equation to be linear but we demonstrate that it is possible to solve non-linear problems using an iterative scheme.

In the special case where multiple measurements are recorded at the same location, i.e. Li

=

L1, i

=

2

,

3

, . . . ,

m, our

method is equivalent to computing the mean valueg

¯

=

(g1

+

g2

+ · · · +

gm)

/

m, and treatingg as a single measurement. By

¯

taking the mean value we obtain a single measurement with a smaller random error. In the case where the measurement locations differ the method can be seen as a way of making the best use of each individual measurement.

Numerical experiments show that the method works well. In particular we apply the method to an application with experimental data [3]. Temperatures at several locations inside a steel slab have been recorded during a heating experiment. In this case the problem is non-linear but with an iterative procedure we are able to solve the problem and compute both the surface temperature and the surface heat-flux.

Our tests also show that a simple method, such as the L-curve, works well for finding the appropriate level of regularization. For our application the rate of decay of the singular values

{

σ

_i

}

, for the operator K1, is known [12] and

thus a formal proof that the proposed method is a regularization of the original problem could be carried out using, e.g., the discrepancy rule as the parameter selection strategy, and possibly an error estimate could be obtained.

There are several potential improvements to the method that could be considered. For instance, the choice of the penalty term in the Tikhonov regularization has a large impact on the result. For each application the appropriate choice must be made. Alternatives include norms

λ

2

_∥

_f

_∥

2

2, first derivatives

λ

2

∥

f

′

_∥

2

2, and second derivatives

λ

2

∥

f

′′

(t)

∥

2 2.

Further, if the primary goal is to compute the surface heat-flux, rather than the temperature, it might work better to instead write the method in terms of an operator H1

:

h1(t)

=

κ

Tx(0

,

t)

↦→

g1(t)

=

T (L1

,

t). This would allow us to directly

apply a penalty term to force smoothness of the computed heat-flux. This is something we will consider in the future as we work on other more real world applications.

(12)

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

Funds for the publication of this paper were provided by the International Science Programme (ISP), through the Eastern Africa Universities Mathematics Programme (EAUMP), University of Rwanda node. The authors are grateful to Patrik Wikström, SSAB, for providing the experimental data used in this paper. The research of Japhet Niyobuhungiro was supported in part by The World Academy of Sciences (TWAS), RGA No. 17-356 RG/MATHS/AF/AC_I - FR3240297744.

References

[1] Chen Hongchu, Frankel Jay I, Keyhani Majid. Nonlinear inverse heat conduction problem of surface temperature estimation by calibration integral equation method. Numer Heat Transfer B 2018;73(5):263–91.

[2] Myrick Justin, Keyhani Majid, Frankel JI. Determination of surface heat flux and temperature using in-depth temperature data. Experimental verification. Int J Heat Mass Transfer 2017;111:982–98.

[3] Wikström Patrik, Blasiak Wlodzimierz, Berntsson Fredrik. Estimation of the transient surface temperature, heat flux and effective heat transfer coefficient of a slab in an industrial reheating furnace by using an inverse method. Steel Res Int 2007;78(1):63–70.

[4] Berntsson Fredrik. An inverse heat conduction problem and improving shielded thermocouple accuracy. Numer Heat Transfer A 2012;61(10):754–63.

[5] Aamir Muhammad, Liao Qiang, Zhu Xun, Rehman Aqeel, Wang Hong, Zubair Muhammad. Estimation of surface heat flux and surface temperature during inverse heat conduction under varying spray parameters and sample initial temperature. Sci World J 2014;2014:721620, p. 13. [6] Eldén Lars, Berntsson Fredrik, Reginska Teresa. Wavelet and fourier methods for solving the sideways heat equation. SIAM J Sci Comput

2000;21(6):2187–205.

[7] Mejia Carlos, Murio DA. Numerical solution of generalized IHCP by discrete mollification. Comput Math Appl 1996;32:33–50. [8] Duda Jan Taler Piotr. Solving direct and inverse heat conduction problems. Springer; 2006.

[9] Berntsson Fredrik. A spectral method for solving the sideways heat equation. Inverse Problems 1999;15(4):891–906.

[10] Cui Miao, Duan Wei-wei, Gao Xiao-wei. A new inverse analysis method based on a relaxation factor optimization technique for solving transient nonlinear inverse heat conduction problems. Int J Heat Mass Transfer 2015;90:491–8.

[11] Lamm Patricia K. Approximation of ill-posed Volterra problems via predictor-corrector regularization methods. Soc Ind Appl Math 1996;56(2):524–41.

[12] Lamm Patricia K, Elden Lars. Numerical solution of frist-kind Volterra equations by sequential tikhonov regularization. Soc Ind Appl Math 1997;34(4):1432–50.

[13] Burggraf O. An exact solution of the inverse problem in heat conduction theory and applications. J Heat Transfer 1964;86(3):373–80. [14] Gustafsson B, Kreiss H-O, Oliger J. Time dependent problems and difference methods. New York: Wiley Interscience; 1995. [15] Hansen Per C. Rank-deficient and discrete ill-posed problems. Soc Ind Appl Math; 1998.

[16] Engel Heinz W, Hanake M, Neubaur A. Regularization of inverse problems. Kruwer Academic Publishers; 1996.

[17] Aster Richard C, Borchers Brian, Thurber Clifford H, editors. 5 Tikhonov regularization. In: Parameter estimation and inverse problems. International geophysics, vol. 90, Academic Press; 2005, p. 89–118.