Simultaneous estimation of spatially distributed thermal conductivity, heat capacity and surface heat transfer coefficient in thermal tomography

International Journal of Heat and Mass Transfer 55 (2012) 7958–7968

Contents lists available at SciVerse ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Simultaneous estimation of spatially distributed thermal conductivity, heat capacity and surface heat transfer coefficient in thermal tomography J.M. Toivanen a, V. Kolehmainen a,⇑, T. Tarvainen a,b, H.R.B. Orlande c, J.P. Kaipio a,d a

Department of Applied Physics, University of Eastern Finland, P.O. Box 1627, 70211 Kuopio, Finland Department of Computer Science, University College London, Gower Street, London WC1E 6BT, UK c Department of Mechanical Engineering, Politecnica/COPPE, Federal Univesity of Rio de Janeiro, Cid. Universitaria, Cx. Postal 68503, 21941-972 Rio de Janeiro, Brazil d Department of Mathematics, University of Auckland, Private Bag 92019, Auckland Mail Centre, Auckland 1142, New Zealand b

a r t i c l e

i n f o

Article history: Received 29 March 2012 Received in revised form 8 August 2012 Accepted 12 August 2012 Available online 10 September 2012 Keywords: Thermal tomography Bayesian inversion Thermal conductivity Heat capacity Surface heat transfer coefficient Non-destructive testing

a b s t r a c t In thermal tomography, the thermal properties of a target are estimated as spatially distributed parameters based on non-invasive measurements of surface temperatures. In the measurement setup, the target is sequentially heated at different source locations and the induced temperature evolutions are measured at several measurement locations on the surface. In [V. Kolehmainen, J. Kaipio, H. Orlande, Reconstruction of thermal conductivity and heat capacity using a tomographic approach, Int. J. Heat Mass Transfer 50 (25–26) (2007) 5150–5160], it was demonstrated with simulations that simultaneous estimation of spatially distributed thermal conductivity and volumetric heat capacity from transient boundary data is feasible when the boundary heat flux from the target to the surrounding medium is known all over the target boundary. In this article, we extend the computational methods towards the more practical setup of imaging targets, where the boundary heat flux from the target to the surrounding medium is not known. We model the surface heat transfer coefficient as a spatially distributed parameter on the target boundary and estimate it simultaneously with the spatially distributed thermal conductivity and volumetric heat capacity using the statistical (Bayesian) inversion framework. The feasibility of the approach is evaluated with simulations. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction In thermal tomography, the unknown thermal conductivity and volumetric heat capacity of a target are estimated as spatially distributed parameters based on temperature or heat flux measurements from the surface of the target. The estimated distributions can be used to infer the composition and structural integrity of the target under inspection. Therefore, thermal tomography can be used for non-destructive testing of targets that withstand moderate temperature changes. The testing can potentially reveal cracks, voids, delamination faults, air bubbles and other manufacturing defects. In the measurement process, the target is heated at a source location and the induced temperature evolutions are measured at multiple measurement locations on the surface. The same heating and measurement process is repeated for a number of source locations. The heatings can be carried out with contact heaters, such as heating resistors that are attached to the surface of the target, and the temperature measurements with contact sensors, such as thermocouples or thermistors. Alternatively, non-contact methods ⇑ Corresponding author. Tel.: +358 40 355 2054. E-mail address: Ville.Kolehmainen@uef.fi (V. Kolehmainen). 0017-9310/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2012.08.024

such as inductive or laser heatings and infrared thermography measurements can be used. The estimation of spatially distributed thermal conductivity and volumetric heat capacity from surface temperature measurements is an ill-posed inverse problem which is unstable with respect to measurement and modeling errors, and requires specific solving methods. In this paper, the inverse problem is considered in the statistical (Bayesian) inversion framework [1,2]. The feasibility of thermal tomography has this far been studied with simulations only. A common feature in all the previous studies is that the surface heat flux from the target to the surrounding medium is assumed known, implying that the boundary condition of the heat equation can be treated as known all over the target boundary. Two-dimensional thermal conductivity distributions were estimated from simulated boundary measurements of the steady state heat transfer in [3]. Transient heat transfer measurements were used in [4,5] and either the thermal conductivity or the heat capacity of the two-dimensional target was estimated, while the other coefficient was assumed to be known. In these works, the target was assumed thermally insulated from the surrounding medium, implying that the surface heat flux can be approximated zero outside the heated part of the boundary. In contrast with [3–5], in [6]

7959

J.M. Toivanen et al. / International Journal of Heat and Mass Transfer 55 (2012) 7958–7968

Nomenclature

X @X mN

Nk qh(x, t) theat mn nj tmeas T(x, t) j(x) c(x) ^ n q(x, t) x t T0

target domain boundary of the domain X number of heater elements location of heater element k heat flux from a heater [W/m2] duration of a heating [s] number of temperature sensors location of temperature sensor j duration of temperature measurements [s] temperature [K] thermal conductivity [W/(m K)] volumetric heat capacity [J/(m3 K)] outward pointing unit normal of @ X heat flux [W/m2] position vector in X [m] time [s] initial temperature of the domain X [K]

the simulated target was not assumed to be thermally insulated from its surroundings. However, the surface heat transfer coefficient determining the surface flux to the surrounding medium was treated as a known parameter while both the two-dimensional thermal conductivity and heat capacity distributions were estimated. A three-dimensional version of thermal tomography, where the thermal conductivity is estimated while the heat capacity and surface heat transfer coefficient are treated as known parameters, was presented in [7]. In [8], the surface heat flux from the target to the surrounding medium was treated as known, and time dependent two-dimensional thermal conductivity distributions were estimated. The shape and thermal conductivity of an inclusion were identified in [9] while the thermal conductivity of the non-inclusion part was assumed to be known. The problem was considered under steady state conditions and either the heat flux or the temperature was known on all of the boundary. Two-dimensional thermal conductivity distributions of a thermal interface were estimated in [10] using multiple different image reconstruction algorithms in such a way that the surface heat transfer coefficient was assumed to be known. In [11], thermal infrared tomography was reviewed and in [12,13] three-dimensional thermal effusivity distributions were obtained by one-sided pulsed thermal imaging. In [11–13], the surface heat transfer coefficient was not considered. In [14], thermal conductivity and volumetric heat capacity were simultaneously estimated as distributed parameters from transient measurement data using the statistical inversion framework. The simulated measurement data was computed from a two-dimensional target which is thermally insulated from its surroundings and thus the surface heat flux was modeled as zero in the non-heated part of the surface. The construction of a practical measurement setup where the surface heat flux is known all over the boundary can be in many applications tedious, or even impossible. Hence, in this paper we extend the approach in [14] to the more practical case where the surface heat flux from the target to the surrounding medium is not known. We model the surface heat transfer coefficient in the boundary condition as a spatially distributed parameter at the boundary and estimate it simultaneously with the spatially distributed thermal conductivity and volumetric heat capacity. The feasibility of the approach is evaluated with simulated measurement data. The rest of the paper is organized as follows. The forward problem of thermal tomography is discussed in Section 2. The inverse problem of thermal tomography is discussed in Section 3. The mea-

h(x)

v(x) Nn

uk(x) ak(t) Nb Dt F ðÞ y

j⁄ c⁄ h⁄ kp cp hp qp e

surface heat transfer coefficient [W/(m2 K)] test function number of nodes in the FE discretization basis functions k in the FE approximations coefficients of the FE approximation of T(x, t) number of boundary nodes in the FE discretization time step of the implicit Euler iteration [s]. FEM based forward problem solution temperature measurements [K] expected values of thermal conductivity j expected values of volumetric heat capacity c expected values of surface heat transfer coefficient h values of k at node p value of c at node p in the FE approximation value of h at node p in the FE approximation value of q at node p in the FE approximation measurement errors [K]

surement setup in the simulations and estimates that are computed are discussed in Section 4 and the results are given in Section 5. Section 6 gives the conclusions. 2. Forward problem 2.1. Forward model Let X  Rn ; n ¼ 2; 3, model the target domain with boundary @ X, Nk  @ X (k = 1, . . . , mN) model the surface patches that are covered by the heater elements and nj 2 @ X(j = 1, . . . , mn) denote the locations of the pointlike temperature sensors. In the measurement process, one of the heaters is turn on at a time and produces a known heat flux qh into the target at Nk during the period t 2 [0, theat]. The resulting temperatures are measured every Dtm seconds during the time interval t 2 [0, tmeas] at all of the measurement locations. This process is then repeated for all mN source locations. Heat transfer is modeled with the initial boundary value problem

@Tðx; tÞ ¼ r  ðjðxÞrTðx; tÞÞ; @t @Tðx; tÞ jðxÞ ^ ¼ qðx; tÞ; x 2 @ X @n Tðx; 0Þ ¼ T 0 cðxÞ

x2X

ð1Þ ð2Þ ð3Þ

where c(x) is the volumetric heat capacity, j(x) is the thermal conductivity, T(x, t) is the temperature, q(x, t) is the heat flux, x is the po^ is the outward pointing unit normal of sition vector in X, t is time, n @ X and T0 is the initial temperature of the target. The heat flux at the surface of the target is written as

qðx; tÞ ¼ qh ðx; tÞ  qc ðx; tÞ;

x 2 @X

ð4Þ

where qh(x, t) is the heat flux imposed on the target by the heater element and qc(x, t) represents the heat loss from the target to the surrounding medium. When heating at Nk, we have

qh ðx; tÞ ¼



qin ; x 2 Nk ; 0;

t 2 ½0; t heat 

otherwise

ð5Þ

where qin is the heat flux from the heater and

qc ðx; tÞ ¼ hðxÞðTðx; tÞ  T 1 Þ;

x 2 @ X; t > 0

ð6Þ

where h(x) is the spatially distributed surface heat transfer coefficient at @ X and T1 is the temperature of the surrounding medium.

7960

J.M. Toivanen et al. / International Journal of Heat and Mass Transfer 55 (2012) 7958–7968

2.2. Finite element approximation of the forward model In this paper, the finite element method (FEM) is used for numerical approximation of the heat transfer problem (1)–(3). We follow the approach described in [14], with some modifications necessitated by the addition of the convective heat flux (6) into the modeling. In the construction of the FE approximation, the heat transfer problem (1)–(3) is written in the variational form

Z @Tðx; tÞ cðxÞ v ðxÞdx ¼ ðqh ðx; tÞ  hðxÞðTðx; tÞ  T 1 ÞÞv ðxÞdS @t X @X Z ð7Þ  jðxÞrTðx; tÞ  rv ðxÞdx

Z

X

Tðx; 0Þ ¼ T 0

ð8Þ

where the test functions v(x) 2 H1(X) and H1(X) is the associated Sobolev space. The domain X is discretized into Ne triangular elements connected at Nn vertex nodes. Using this mesh, the solution of the variational form (7) is approximated as

Tðx; tÞ 

Nn X

ak ðtÞuk ðxÞ 2 Q h

ð9Þ

k¼1

Nn X

jp up ðxÞ

ð10Þ

p¼1 Nn X cp up ðxÞ

cðxÞ 

ð11Þ

p¼1

X

hðxÞ 

hp up ðxÞj@ X

ð12Þ

p2@ X

qh ðx; tÞ 

X

qp ðtÞup ðxÞj@ X

ð13Þ

p2Nk

where jp, cp, hp and qp(t) are the values of thermal conductivity, volumetric heat capacity, surface heat transfer coefficient and heat flux at node p, and p 2 @ X are the indices of basis functions on @ X and p 2 Nk on the active heater. Inserting the approximations (9)–(13) into (7) and choosing v(x) = uj(x), we get the matrix form

M

@a ¼ Ga  F a þ V þ Q @t

Z X Nn cp up ðxÞuk ðxÞuj ðxÞdx

Vj ¼

Z X Nn Z

Z

jp up ðxÞruk ðxÞ  ruj ðxÞdx

X

@ X p2N k

Qj ¼

Z

ð16Þ

hp up ðxÞuk ðxÞuj ðxÞdS

ð17Þ

qp ðtÞup ðxÞuj ðxÞdS

ð18Þ

@ X p2@ X

X

Given the vectors al, approximations for the temperatures T(x, t) at the measurement locations at measurement instants are obtained from (9) and interpolation in the time domain. In the following, we use the notation

F ðj; c; hÞ for the FEM based forward problem solution, where j 2 RNn ; c 2 RNn and h 2 RNb are vectors of coefficients of j(x), c(x) and h(x) in Eqs. (10)–(12) and Nb is the number of boundary nodes. The measured data in the thermal tomography experiment is arranged in measurement vector y 2 Rmt mN mn as

 T y ¼ T ð1Þ ; T ð2Þ ; . . . ; T ðmN Þ

ð21Þ

X

@ X p2@ X

hp up ðxÞT 1 uj ðxÞdS

where j, k = 1, . . . , Nn.

where T 2 Rmt mn are the temperature recordings corresponding to heating k so that

T ðkÞ ¼ ðTðn1 ; t1 Þ; . . . ; Tðn1 ; t meas Þ; Tðn2 ; t 1 Þ; . . . ; Tðnmn ; t meas ÞÞT

ð22Þ

where mt is the number of measurement times. The random measurement errors are modeled as additive zero mean Gaussian errors which are mutually independent with the unknowns. Thus, we write the measurement model for the thermal tomography experiment as

y ¼ F ðj; c; hÞ þ e;

e  N ð0; Ce Þ

ð23Þ

where Ce is the covariance matrix of the measurement errors e. 3. Inverse problem 3.1. Posterior density model

pðj; c; hjyÞ ¼

X p¼1

F jk ¼

al ¼ ½a1 ðlDtÞ; . . . ; aNn ðlDtÞT :

ð15Þ

and the elements of the matrices and vectors are

X p¼1

where Dt is the time step and

ð14Þ

a :¼ aðtÞ ¼ ½a1 ðtÞ; . . . ; aNn ðtÞT

Gjk ¼

ð20Þ

In this paper, the thermal tomography inverse problem is the simultaneous estimation of spatially distributed thermal conductivity, volumetric heat capacity and surface heat transfer coefficient from boundary temperature measurements. The inverse problem is solved using the statistical (Bayesian) inversion framework [1,2]. In statistical inversion, all unknown and measured parameters are modeled as random variables and the uncertainty related to their values is encoded into their probability density models. The solution of the inverse problem is the posterior probability density p(j, c, hjy), given by the Bayes’ theorem as

where

M jk ¼

alþ1 ¼ ðI þ M1 ðG þ FÞDtÞ1 ðal þ M1 ðV þ Q ÞDtÞ

ðkÞ

where ak(t) are the coefficients of the approximation, uk(x) are the piecewise linear nodal basis functions and the space Qh = span{uk}  H1(X). For the discretization of the thermal parameters, we use the same finite dimensional basis as for the temperatures, leading to

jðxÞ 

For the solution of the semidiscrete form (14), we employ the implicit Euler method with constant time step. This gives the solution at time (l + 1)Dt by

ð19Þ

pðyjj; c; hÞpðj; c; hÞ pðyÞ

ð24Þ

where p(yjj, c, h) is the likelihood model, p(j, c, h) is the prior model and p(y) is the marginal density of the measurement, which plays the role of a normalization constant. The posterior model represents the complete probabilistic model of the inverse problem given the measured data and a priori information that is encoded in the prior model. The likelihood model is a statistical model for the measurement and it describes the probability of different measurement realizations for a given realization of (j, c, h). Given the measurement model (23) and the assumptions on e, the likelihood can be written as

7961

J.M. Toivanen et al. / International Journal of Heat and Mass Transfer 55 (2012) 7958–7968

2

pðyjj; c; hÞ ¼ pe ðy  F ðj; c; hÞÞ

  1 / exp  kLe ðy  F ðj; c; hÞÞk2 2

ð25Þ

where pe is the probability density of the measurement errors e, and LTe Le ¼ C1 e , where Ce is the covariance matrix of the measurement errors. The prior model is a statistical model for the a priori information. In this study, we model (j, c, h) as mutually independent Gaussian random variables. Thus, for each parameter we use a marginal prior model of the form

  1 pðf Þ / exp  kLf ðf  f Þk2 2

ð26Þ

where f is either j, c or h, f⁄ is the corresponding expected value, Cf the covariance matrix and LTf Lf ¼ C1 f . Thus the joint prior model becomes

pðj; c; hÞ ¼ pðjÞpðcÞpðhÞ

 1 / exp  ðkLj ðj  j Þk2 þ kLc ðc  c Þk2 þ kLh ðh  h Þk2 Þ : 2

Lc

0

0

By Eqs. (24), (25) and (27), the posterior model is

 1 pðj; c; hjyÞ / exp  kLe ðy  F ðj; c; hÞÞk2 þ kLj ðj  j Þk2 2    2 ð28Þ þ Lc ðc  c Þk þ kLh ðh  h Þk2 :

In statistical inversion, the solution of the inverse problem is analyzed and visualized by point and spread estimates from the posterior model. Two most common point estimates are the conditional mean (CM) and the maximum a posteriori (MAP) estimate [1]. The computation of the CM estimate and the spread estimates is an integration problem, which in large dimensional non-linear problems usually requires the use of Markov chain Monte Carlo (MCMC) integration methods [2,15,16]. The MCMC methods are based on random sampling of the posterior density and the sampling process requires a very large number of forward solutions. For problems such as thermal tomography where the forward solution involves solving FEM approximation of a PDE, they lead to very long computation times. The computation of the MAP estimate is an optimization problem and as such computationally significantly less expensive than the integration based estimates. In this study, we compute the MAP estimate which is obtained by

n ðj; c; hÞMAP ¼ arg min kLe ðy  F ðj; c; hÞÞk2 þ kLj ðj  j Þk2 j;c;h o þkLc ðc  c Þk2 þ kLh ðh  h Þk2 :

ð32Þ

J ¼ ½ Jj

Jc

Jh 

ð33Þ

where the blocks are

2 6 6 Jj ¼ 6 4 2 6 6 Jc ¼ 6 4 and

2

6 6 Jh ¼ 6 6 4

@T ð1Þ @ j1

...

3

@T ð1Þ @ jN n

.. .

.. .

ðmN Þ

ðmN Þ

@T @ j1

...

@T @ jN n

@T ð1Þ @c1

...

@T ð1Þ @cNn

.. .

@T ðmN Þ @c1

...

@T ðmN Þ @cNn

@T ð1Þ @h1

...

@T ð1Þ @hN b

.. .

@T ðmN Þ @h1

@T ðmN @hN

...

ðj; c; hÞðiþ1Þ ¼ ðj; c; hÞðiÞ þ ^sðiÞ dðj; c; hÞðiÞ

7 7 7 5

2R

is obtained as the

  h   J T LTe Le J þ LT L dðj;c;hÞðiÞ ¼ J T LTe Le y  F ððj;c;hÞði1Þ Þ  T  ði1Þ ð31Þ þLT L j  jði1Þ ;c  cði1Þ ;h  h where J :¼ J((j, c, h)(i)) is the Jacobian matrix,

ð35Þ

3 7 7 7: 7 5 Þ

ð36Þ

b

The elements of the blocks of the Jacobian can be derived from the implicit Euler iteration (20) by taking the appropriate partial derivatives. This gives iterations

@ alþ1 @ al @M 1 @M1 ¼ B1  ðF þ GÞDt alþ1 þ ðV þ Q ÞDt @ck @ck @ck @ck  @ alþ1 @ al @Q @F ¼ B1 þ M 1 Dt M 1 Dt alþ1 @hk @hk @hk @hk

ð37Þ ! ð38Þ ð39Þ

where

B1 ¼ ðI þ M 1 ðF þ GÞDtÞ1

ð40Þ

@M 1 @M 1 ¼ M1 M : @ck @ck

ð41Þ

In the same way as for the temperatures in (20), interpolation can be done in the time domain to get the partial derivatives at measurement instants. 4. Methods

In all the simulation cases in this paper, the target domain

ð29Þ

ð30Þ 2N n þN b

ð34Þ

3

.. .

.. .

7 7 7 5

4.1. Simulation of the measurements & test cases

We compute the MAP solution (29) with Gauss–Newton optimization, where the estimate is iteratively updated as

where the search direction dðj; c; hÞ solution of

7 05 Lh

 @ alþ1 @ al @G ¼ B1  M1 Dt alþ1 @ jk @ jk @ jk

3.2. Computation of the MAP-estimate

ðiÞ

0

and the step lengths ^sðiÞ are determined by a line search algorithm [17]. The Jacobian is of the form



ð27Þ

0

Lj 6 L¼4 0

3

X  R2 is a circle with radius of 2 cm. The target is heated from its boundary with heater elements at patches Nk, k = 1, . . . , 8 and each heater element covers 1/8 of the boundary so that [8k¼1 Nk ¼ @ X.

The heaters produce one at a time a heat flux qin = 1000 W/m2 for theat = 240 s and temperature evolutions are recorded by pointlike temperature sensors at locations nj, j = 1, . . . , 8. The temperature sensors are located between the heater elements and record temperatures for tmeas = 2000 s every Dtm = 4 s starting from the beginning of a heating. Thus, the measurement vector y 2 R32;000 . The targets are surrounded by air at temperature T1 = 23 °C which is also the initial temperature T0 of the target. An illustration of the measurement geometry is shown in Fig. 1. We consider two test cases. The first test case (‘‘case I’’) simulates a simple case of non-destructive testing from a homogeneous

7962

J.M. Toivanen et al. / International Journal of Heat and Mass Transfer 55 (2012) 7958–7968

material. Here the thermal properties correspond to a homogeneous clay object containing a circular air cavity as an imperfection. The target distributions (j, c) are shown in the top row in Fig. 3. The surface heat transfer coefficient h corresponds to a spatially distributed coefficient caused by buoyancy induced flow [18,19] and is shown with a dashed line in Fig. 4. The second test case concerns imaging of a non-homogeneous target structure. This test case is divided into two parts, namely ‘‘case IIa’’ and ‘‘case IIb’’. In case IIa, the domain consist of an unknown layered structure with thermal properties corresponding to two types of clay. The target distributions are shown in the first and second columns on the first row in Fig. 9. In case IIb, the domain is otherwise as in case IIa but contains an air cavity as an imperfection to illustrate non-destructive testing from a nonhomogeneous target structure. The target distributions for case IIb are shown in the third and fourth columns on the first row in Fig. 9. The surface heat transfer coefficient is in cases IIa and IIb the same piecewise constant function where the value of h depends locally on the type of the clay. The surface heat transfer coefficient is shown with a dashed line in Fig. 10. Thermal properties used in cases I, IIa and IIb are shown in Table 1. The discretization details of the finite element meshes that were used for computation of the simulated measurement data in cases I, IIa, IIb and in all of the inverse problems are given in Table 2. In the selection of the finite element mesh for the inverse problem, we studied the convergence of the FEM approximation by computing solutions in a sequence of meshes with fixed (j, c, h) and compared the solutions against a reference solution computed in a very accurate discretization with Nn = 10127 node points. The solution was found practically converged at 500 nodes, and the level of discretization error at this discretization level was negligible compared to the level of the measurement noise. We selected a mesh with Nn = 581 nodes for the computations. To avoid the so-called inverse crime, the simulated measurement data was computed using finer meshes than are used in the estimation. The noisy realizations of the simulated measurement data were obtained by adding Gaussian zero mean random noise with standard deviation re = 0.1 °C to the simulated noiseless measurements. Example of simulated measurement data is shown in Fig. 2, which shows the measured temperature evolutions for a temperature sensor next and opposite to the active heater in the first heating in case I.

Table 1 Values of thermal conductivity j (W/(m K)), volumetric heat capacity c (J/(m3 K)) and surface heat transfer coefficient h (W/(m2 K)) for the test cases. Case

Material 1

I IIa IIb ⁄

Material 2

c

h

j

c

h

j

c

1 1 1

106 106 106

⁄

– 1.25 1.25

– 1.25  106 1.25  106

– 2 2

0.03 – 0.03

1000 – 1000

10 10

See Fig. 4.

Table 2 The numbers of elements Ne, nodes Nn and boundary nodes Nb in the triangular FEM meshes for the simulation of measurements in cases I, IIa, IIb and solving of the inverse problems. Case

Ne

Nn

Nb

I IIa/b Inverse

2400 1938 1096

1269 1036 581

136 132 64

4.2. Initial estimate for the MAP estimation To ensure convergence and obtain the MAP estimate in a reasonably small number of Gauss–Newton iterations (30), the initial estimate for the iteration should be carefully selected. A widely used procedure for finding the initial estimate in other diffusion equation based tomography modalities such as electrical impedance tomography and diffuse optical tomography is to estimate the average values of the distributed parameters that give best match to the measured data in the least squares sense. Within the statistical inversion framework, this initial estimation procedure corresponds to parameterizing the unknowns as constant functions

jðxÞ ¼

Nn X

j0 up ðxÞ

ð42Þ

p¼1 Nn X c0 up ðxÞ

cðxÞ ¼

ð43Þ

p¼1

X

hðxÞ ¼

h0 up ðxÞj@ X

ð44Þ

p2@ X

measurement sensor 1

35

Ξ

3π/4

3

ξ

ξ2

4

20

ξ

1

ξ

6

−3π/4

Ξ

Ξ6

ξ ξ7

7

200

400

600

800

1000 1200 1400 1600 1800 2000

Ξ8 −π/4

25 24 23 22

0

−π/2 Fig. 1. The measurement geometry showing the location of the heaters at patches Nk, k = 1, . . . , 8 and the temperature sensors at location nj, j = 1, . . . , 8. Also the surface angles h (rad) used in the visualization of the surface heat transfer coefficients are shown.

measurement sensor 5

26

0

8

Ξ

0

t [s]

T [ oC]

5

25

1

5

Ξ

π/4

2

ξ

±π

30

o

Ξ

3

ξ

Ξ4

T [ C]

π/2

Inclusion

j

200

400

600

800

1000 1200 1400 1600 1800 2000

t [s] Fig. 2. Example measurement data from case I for the first heating. Temperatures T (°C) as a function of time t (s). Top: Data from a sensor next to the active heater element. Bottom: Data from a sensor opposite to the active heater element. The dashed line marks the time instant t = 240 s when the heater element is turned off.

7963

J.M. Toivanen et al. / International Journal of Heat and Mass Transfer 55 (2012) 7958–7968

and estimate the values of the scalar coefficients j0 2 R; c0 2 R; h0 2 R by maximizing the likelihood p(yjj0, c0, h0). This estimate is obtained by

n o ^ 0 Þ ¼ arg min kLe ðy  F ðj0 ; c0 ; h0 ÞÞk2 : ^ 0 ; ^c0 ; h ðj

ð45Þ

j0 ;c0 ;h0

We employ the Gauss–Newton method for the solution of (45). Once the problem (45) has been solved, the initial estimate for the MAP estimation is set to constant vectors with values defined ^0 Þ. ^ 0 ; ^c0 ; h by the estimated scalars ðj 4.3. Reference estimates To assess the effect of having the surface heat transfer coefficient h as an additional distributed unknown, we show reference estimates in which j and c are estimated with a fixed realization of h. These estimates are obtained as

n ðj; cÞREF ¼ arg min kLe ðy  F ðj; c; hfixed ÞÞk2 þ kLj ðj  j Þk2 j;c o þkLc ðc  c Þk2

ð46Þ

where hfixed is the fixed realization of h. We consider the following two cases: 1. The surface heat transfer coefficient is known exactly hfixed:¼htrue(x). This represents reference of estimating (j, c) in the ideal setup with an exactly known boundary condition and is in that sense equivalent to the previously published thermal tomography results (see e.g. [14]). 2. Using incorrect realization hfixed :¼ htrue ðxÞ which is chosen as the mean of the true surface heat transfer coefficient. This represents the case of estimating (j, c) with a misspecified boundary condition. Fig. 3. Test case I. Top row: True thermal conductivity j (W/(m K)) (left) and volumetric heat capacity c (J/(m3 K)) (right). Second row: Estimates using the correct fixed surface heat transfer coefficient hfixed = htrue. Third row: Estimates using the incorrect fixed surface heat transfer coefficient hfixed ¼ htrue . Bottom row: Estimates with the simultaneous estimation of a spatially distributed surface heat transfer coefficient.

The initial estimate for the iterative solution of (46) can be computed similarly as above, that is, by (45) with the modification that h is treated as fixed and known. 4.4. Prior models As the prior density model of thermal conductivity and volumetric heat capacity, we use a proper smoothness prior model constructed similarly as in [20,21]. We model f = j,c in two mutually independent parts

11 10

f ¼ fin þ fbg

ð47Þ

9

h [W/(m2K)]

8

Table 3 Maximum (dMAX) and root mean squared (dRMS) errors of the estimates of thermal conductivity j and volumetric heat capacity c in cases I, IIa and IIb. Column with htrue corresponds to errors of parameters estimated using the fixed true value of the surface heat transfer coefficient, htrue using the fixed mean value and est. h(x) to estimates where the surface heat transfer coefficient was simultaneously estimated.

7 6 5 4

I, I, I, I,

3 2 1 −4

true mean distributed

−3

−2

−1

0

1

2

3

4

θ [rad] Fig. 4. Test case I. True (dashed line) and spatially distributed estimate (dot-line) of the surface heat transfer coefficient h (W/m2 K) as a function of surface angle h (rad). The solid line shows the mean htrue of the true surface heat transfer coefficient. Corresponding angles on the surface of the target are shown in Fig. 1.

dRMS(j) dRMS(c) dMAX(j) dMAX(c)

htrue

htrue

Est. h(x)

0.125 1.57  105 0.503 6.07  105

0.316 3.06  105 0.999 10.2  105

0.108 1.78  105 0.458 6.21  105

IIa, IIa, IIa, IIa,

dRMS(j) dRMS(c) dMAX(j) dMAX(c)

0.071 0.71  105 0.200 2.38  105

0.188 1.05  105 0.630 2.71  105

0.080 0.69  105 0.240 2.24  105

IIb, IIb, IIb, IIb,

dRMS(j) dRMS(c) dMAX(j) dMAX(c)

0.104 1.37  105 0.427 0.71  105

0.251 1.98  105 0.948 1.01  105

0.128 1.31  105 0.477 0.70  105

7964

J.M. Toivanen et al. / International Journal of Heat and Mass Transfer 55 (2012) 7958–7968

Fig. 5. True values (first row) and estimates of thermal conductivity (columns 1, 3 and 5) and volumetric heat capacity (columns 2, 4 and 6) for more conductive targets. Estimates on the rows are as in Fig. 3. Background thermal conductivity of the target is (a) 2, (b) 5 and (c) 10 W/(m K).

5

0.25

2.8

x 10

RMSE(c)

RMSE(κ)

2.6 0.2

0.15

2.4 2.2 2

0.1 0

1

2 σe

3

4

1.8 0

1

2 σe

3

4

1

2 σe

3

4

10

0.7

9

max error(c)

max error(κ)

5

0.8

0.6 0.5 0.4 0

1

2 σe

3

4

x 10

8 7 6 0

Fig. 6. RMS- and maximum error of the estimated thermal conductivity and heat capacity in case I as a function of the standard deviation of the measurement noise. Vertical dotted lines indicate the standard deviations 0.67, 2.33 and 4 °C of the noise level in the measurements in Fig. 7 and the three estimates in Fig. 8.

J.M. Toivanen et al. / International Journal of Heat and Mass Transfer 55 (2012) 7958–7968

where fin is zero mean and spatially inhomogeneous, and fbg is spatially constant with non-zero mean. For fbg, we write fbg ¼ gI, where I 2 RNn is a vector of ones and g is a scalar random variable with Gaussian distribution g  N ðf ; r2bg;f Þ. For fin, we write fin  N ð0; Cin;f Þ. Since we model fin and fbg as mutually independent, we have covariance

Cf ¼ Cin;f þ r2bg;f IIT :

ð48Þ

In the construction of the covariance Cin,f of the inhomogeneous part, the approximate correlation lengths are adjusted to match the size of the expected inhomogeneities and the marginal variances of the elements of f are equalized to the given value r2in;f . The correlation length was set as 1 cm for both the thermal conductivity and the volumetric heat capacity. The mean of the prior was set to f ¼ ftrue where ftrue is the mean of the true f = j, c. The standard deviations were set to rin,f = 0.3f⁄ and rbg,f = 0.1f⁄, leading to marginal standard deviation rf = rin,f + rbg,f = 0.4f⁄. With these settings, the three standard deviation (i.e., 99.7%) confidence limits of the prior model are f⁄ ± 1.2f⁄. This can be interpreted such that in the prior density model the values of f are expected to lie within this interval with probability of 99.7%. The values of the surface heat transfer coefficient h are modeled as mutually independent and identically distributed Gaussian random variables. Therefore, we have Ch ¼ r2h Ib , where rh is the standard deviation of the surface heat transfer coefficient and Ib is a Nb Nb identity matrix. The mean and the standard deviation were set as h ¼ htrue and rh = 0.3h⁄, where htrue is the mean of the true surface heat transfer coefficient. Thus, in this prior density model the values of h are expected to lie within h⁄ ± 0.9h⁄ with probability of 99.7%.

sensor 1

7965

5. Results and discussion 5.1. Case I The results for case I with the thermal parameters corresponding to the buoyancy induced flow and the air cavity in the otherwise homogeneous clay domain are shown in Figs. 3 and 4. The top row in Fig. 3 shows the true thermal conductivity (left) and volumetric heat capacity (right). The second row shows the reference estimate (46) with the correct fixed realization hfixed = htrue. The true surface heat transfer coefficient htrue is shown with a dashed line in Fig. 4. The third row shows the estimate (46) with the incorrect fixed realization hfixed ¼ htrue , which is shown with a solid line in Fig. 4. The bottom row shows the estimate (29) in which the surface heat transfer coefficient is estimated simultaneously with the thermal conductivity and volumetric heat capacity. The estimated spatially distributed surface heat transfer coefficient is shown with a dot-line in Fig. 4. The estimates that are obtained using the fixed but incorrect surface heat transfer coefficient (third row in Fig. 3) show artefacts near the boundary, especially in the estimate of thermal conductivity. Furthermore, the location and shape of the air cavity is not observable in either estimate. This is caused by the incorrectness of the boundary condition (3) due to the incorrect fixed realization of the surface heat transfer coefficient. The estimates on the bottom row, in which the surface heat transfer coefficient is estimated simultaneously with the thermal conductivity and volumetric heat capacity, are in contrast free of these artefacts and similar to the

sensor 5

o T [ C]

35 30 25 20 15

T [oC]

35 30 25 20 15

o

T [ C]

35 30 25 20 15 0

1000

t [s]

2000

1000

2000

t [s]

Fig. 7. Measurement data from measurement sensors 1 (next to the active heater, left column) and 5 (opposite to the active heater, right column) for heating 1 with three levels of measurement errors. Top to bottom: re = 0.67, 2.33 and 4 °C. The corresponding estimates of (j, c) are shown in Fig. 8.

Fig. 8. Thermal conductivity (left column) and heat capacity (right column). Top to bottom: True values, estimates with noise levels re = 0.67, 2.33 and 4 °C, respectively. The estimate (case I) using noise level re = 0.1 °C is shown in bottom row in Fig. 3.

7966

J.M. Toivanen et al. / International Journal of Heat and Mass Transfer 55 (2012) 7958–7968

reference estimates in the second row which have been estimated using the correct fixed realization hfixed = htrue. The root mean squared (RMS) errors

dRMS ð^f Þ ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z jftrue ðxÞ  ^f ðxÞj2 dx X

and maximum errors

dMAX ð^f Þ ¼ max jftrue ðxÞ  ^f ðxÞj of the estimates of thermal conductivity and volumetric heat capacity for case I are given in Table 3. The proposed method was also tested with data from targets with more conductive material than used in Fig. 3. The background thermal conductivity in cases (a)–(c) in Fig. 5 is two, five and ten fold compared to the test case in Fig. 3. As can be seen, the proposed approach (bottom row in Fig. 5) performs equally well as the reference where the true surface heat transfer coefficient is used (second row) also when imaging more conductive targets.

An interesting observation from Fig. 5 is that the estimation seems to be less error prone for the use of incorrect fixed value of the surface heat transfer coefficient when the conductivity increases (third row in Fig. 5), leading to almost as good estimates as the reference method (second row) and the proposed approach (fourth row) when the conductivity is ten fold to the original. To evaluate the robustness of the approach with respect to the level of the measurement noise, estimates were computed using the target in case I with the same heating parameters but with higher levels of noise in the data. The RMS and maximum errors of these estimates with respect to the noise level are shown in Fig. 6. Realizations of measurement data at two measurement sensors corresponding to noise levels re = 0.67, 2.33 and 4 °C are shown in Fig. 7 and corresponding estimates of (j, c) in Fig. 8. As can be seen, the proposed approach tolerates noisy data well, even the estimate with noise level re = 4 °C (noise standard deviation 40-fold compared to case I in Fig. 3) is indicative of the fault region. Note that in case of targets that tolerate high temperature rises, the relative noise level (SNR) could in practice be improved by

Fig. 9. Test case IIa (left) and IIb (right). True and estimated values of thermal conductivity j (W/(m K)) and volumetric heat capacity c (J/(m3 K)) for the layered targets. The color maps are scaled to the minimum and maximum values of all estimates. Estimates on the rows are as in Fig. 3.

J.M. Toivanen et al. / International Journal of Heat and Mass Transfer 55 (2012) 7958–7968

7967

5.3. Discussion

12 true mean

10

intact cavity

h [W/(m2K)]

8

6

4

2

0 −4

−3

−2

−1

0

1

2

3

4

θ [rad] Fig. 10. Test case IIa and IIb. True (dashed line), the mean htrue (solid line) and estimated spatially distributed surface heat transfer coefficient h (W/m2 K) as a function of surface angle h (rad) for the intact target (IIa) (black dot-line) and the non-intact target with the air cavity (IIb) (gray dot-line). Corresponding angles on the surface of the target are shown in Fig. 1.

applying more powerful heating into the target when using sensors with high uncertainty.

The estimation of thermal conductivity and volumetric heat capacity with a fixed surface heat transfer coefficient leads to infeasible estimates when the heat transfer coefficient is not accurately known. The proposed approach of estimating the surface heat transfer coefficient as a distributed parameter on the surface simultaneously with the thermal conductivity and volumetric heat capacity gives similar estimates as the reference approach in which the surface heat transfer coefficient is exactly known. These results suggest that the proposed method facilitates thermal tomographic imaging in cases where the boundary condition describing the surface heat flux from the target to the surrounding medium is not known, implying a significantly simpler experimental measurement setup (e.g., no need to insulate the target from the surrounding). In the estimates of thermal conductivity and volumetric heat capacity, where visible, the layer interfaces are somewhat blurred. Such blurring of sharp boundaries is typical behavior of diffuse tomography imaging, especially when smoothness promoting prior models are employed because in such prior models, sharp boundaries have a relatively low prior probability. However, the proposed method performed robustly despite using such ‘‘nonoptimal’’ prior models. In practice, the reconstruction of targets that are expected to contain edges with large jumps in the parameter values could be improved by utilizing edge-preserving prior models such as the total variation prior [22,23]. 6. Conclusions

5.2. Cases IIa and IIb The results for cases IIa and IIb corresponding to the layered clay targets are shown in Figs. 9 and 10. The top row in Fig. 9 shows the true thermal conductivity and volumetric heat capacity (columns 1–2 from the left for case IIa and columns 3–4 for case IIb). The second row shows the reference estimates (46) using the correct fixed realization hfixed = htrue. The true surface heat transfer coefficient htrue is shown with a dashed line in Fig. 10. The third row shows the estimates (46) with the incorrect fixed realization  true , which is shown with a solid line in Fig. 10. The bottom hfixed ¼ h row shows the estimates (29) in which the surface heat transfer coefficient is estimated simultaneously with the thermal conductivity and volumetric heat capacity. The estimated surface heat transfer coefficient is shown with a dot-line (black for case IIa and grey for IIb) in Fig. 10. Using the incorrect surface heat transfer coefficient (third row in Fig. 9) leads to highly erroneous estimates of thermal conductivity. Whereas the clay layers are not visible in the estimates of thermal conductivity, they are distinguished to some extent in the estimates of volumetric heat capacity. However, in case IIb the air cavity is not visible in either estimate when using the incorrect fixed realization of h. Furthermore, the estimates of IIb are almost identical to those of case IIa, making non-destructive testing of the target’s structure unreliable. In contrast, the estimates on the bottom row, in which the surface heat transfer coefficient is estimated simultaneously with the thermal conductivity and volumetric heat capacity, are free of artefacts and are similar to the reference estimates in the second row, which have been estimated using the correct fixed realization hfixed = htrue. In these estimates, the air cavity is clearly seen in both estimates. Furthermore, the presence of the cavity has not significantly affected the estimate of the surface heat transfer coefficient which ideally should be the same in both cases. The RMS and maximum errors of the estimates of thermal conductivity and volumetric heat capacity for cases IIa and IIb are given in Table 3.

In this paper, we evaluated the feasibility of a computational approach for thermal tomographic imaging in cases where the boundary condition describing the surface heat flux from the target to the surrounding medium is not known. In the proposed approach, thermal conductivity, volumetric heat capacity and surface heat transfer coefficient are simultaneously estimated as distributed parameters using the statistical inversion framework. To evaluate the approach, we examined simulation cases with spatially distributed surface heat transfer coefficients and different internal structures. The proposed method was found to produce similar estimates of thermal conductivity and volumetric heat capacity as are obtained when the surface heat transfer coefficient is exactly known, implying that the extension of thermal tomography to imaging of targets where the heat flux from the target to the surrounding medium is not known is feasible. Acknowledgements This work was supported by the Academy of Finland (Projects 136220, 140984, 119270, 140731, 213476 and 250215, Finnish Centre of Excellence in Inverse Problems Research) and the Magnus Ehrnrooth foundation. The support for the visits of Ville Kolehmainen and Jari Kaipio to COPPE/UFRJ, provided by CNPq, CAPES and FAPERJ, is also greatly appreciated. References [1] J. Kaipio, E. Somersalo, Statistical and Computational Inverse Problems, Springer Science+Business Media, Inc., 2005. [2] J. Kaipio, C. Fox, The Bayesian framework for inverse problems in heat transfer, Heat Transfer Eng. 32 (9) (2011) 718–753. [3] M. Jones, A. Tezuka, Y. Yamada, Thermal tomographic detection of inhomogeneities, J. Heat Transfer 117 (1995) 969. [4] R. Kline, W. Winfree, V. Bakirov, A new approach to thermal tomography, AIP Conference Proceedings 657 (2003) 682. [5] V. Bakirov, R. Kline, W. Winfree, Discrete variable thermal tomography, AIP Conference Proceedings 700 (2004) 469. [6] V. Bakirov, R. Kline, W. Winfree, Multiparameter thermal tomography, AIP Conf. Proc. 700 (2004) 461.

7968

J.M. Toivanen et al. / International Journal of Heat and Mass Transfer 55 (2012) 7958–7968

[7] V. Bakirov, R. Kline, Recent advances in thermal tomography, AIP Conference Proc. 760 (2005) 875. [8] C. Huang, S. Chin, A two-dimensional inverse problem in imaging the thermal conductivity of a non-homogeneous medium, Int. J. Heat Mass Transfer 43 (22) (2000) 4061–4071. [9] M. Ardakani, M. Khodadad, Identification of thermal conductivity and the shape of an inclusion using the boundary elements method and the particle swarm optimization algorithm, Inverse Probl. Sci. Eng. 17 (7) (2009) 855– 870. [10] H. Erturk, Evaluation of image reconstruction algorithms for non-destructive characterization of thermal interfaces, Int. J. Therm. Sci. 50 (2011) 906–917. [11] V. Vavilov, D. Nesteruk, V. Shiryaev, A. Ivanov, W. Swiderski, Thermal (infrared) tomography: terminology, principal procedures, and application to nondestructive testing of composite materials, Russ. J. Nondest. Test. 46 (3) (2010) 151–161. [12] J. Sun, Quantitative three-dimensional imaging by thermal tomography method, AIP Conf. Proc. 1335 (2011) 430. [13] J. Sun, D. Thompson, D. Chimenti, Quantitative thermal tomography imaging of complex material structures, AIP Conf. Proc. 1430 (1) (2012) 507. [14] V. Kolehmainen, J. Kaipio, H. Orlande, Reconstruction of thermal conductivity and heat capacity using a tomographic approach, Int. J. Heat Mass Transfer 50 (25–26) (2007) 5150–5160.

[15] L. Tierney, Markov chains for exploring posterior distributions, Ann. Stat. (1994) 1701–1728. [16] A. Smith, G. Roberts, Bayesian computation via the Gibbs sampler and related Markov chain Monte Carlo methods, J. Roy. Stat. Soc. Ser. B (Meth.) (1993) 3– 23. [17] J. Nocedal, S.J. Wright, Numerical Optimization, second ed., Springer, New York, 2006. [18] A. Bejan, A. Kraus, Heat Transfer Handbook, vol. 1, Wiley-Interscience, 2003. [19] T. Misumi, K. Suzuki, K. Kitamura, Fluid flow and heat transfer of natural convection around large horizontal cylinders: experiments with air, Heat Transfer–Asian Res. 32 (4) (2003) 293–305. [20] V. Kolehmainen, M. Schweiger, I. Nissilä, T. Tarvainen, S. Arridge, J. Kaipio, Approximation errors and model reduction in three-dimensional diffuse optical tomography, JOSA A 26 (10) (2009) 2257–2268. [21] V. Kolehmainen, T. Tarvainen, S. Arridge, J. Kaipio, Marginalization of uninteresting distributed parameters in inverse problems – application to optical tomography, Int. J. Uncert. Quant. 1 (1) (2011) 1–17. [22] L.I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60 (1992) 259–268. [23] J. Kaipio, V. Kolehmainen, E. Somersalo, M. Vauhkonen, Statistical inversion and Monte Carlo sampling methods in electrical impedance tomography, Inverse Probl. 16 (2000) 1487.