Identification of temperature dependent thermal conductivity using the gradient method

IDENTIFICATION OF TEMPERATURE DEPENDENT THERMAL CONDUCTIVITY USING THE GRADIENT METHOD

Ewa Majchrzak^1,2, Katarzyna Freus², Sebastian Freus²

1 Silesian University of Technology, Gliwice, Poland

2 Institute of Mathematics and Computer Science, Czestochowa University of Technology, Poland

Abstract. In the paper the application of the BEM for numerical solution of the inverse parametric problem is presented. On the basis of the knowledge of temperature field in the domain considered the temperature dependent thermal conductivity is identified. The steady state is considered (2D problem) and mixed boundary conditions are taken into account. The inverse problem is solved using the gradient method. In the final part of the paper the results of computations are shown.

1. Formulation of the problem

The following 2D problem is considered

) ( ) ( n ) ( λ ) ( :

) ( ) ( :

0 )]

( ) ( λ [ :

2 1

x q x T T x

q x

x T x T x

x T T x

b b

=

∇

⋅

−

= Γ

∈

= Γ

∈

=

∇

∇ Ω

∈

(1)

where T is the temperature, x = (x1,x2) are the spatial coordinates, λ(T) is the thermal conductivity, Tb(x), qb(x) are known boundary temperature and boundary heat flux, n is the normal outward vector at the boundary point x. We assume that

2

) 1

(

λ T =cT+c (2)

where c1, c2 are the coefficients.

If the direct problem is considered then all geometrical and thermophysical parameters appearing in the mathematical model (1) are known.

In the paper the inverse parametric problem is considered in which we assume that the coefficients c1, c2 are unknown. In order to solve the inverse problem the additional information is necessary. So, we assume that the temperatures at the selected points xi ∈Ω are given

M i

x x T

T_{d i}= _d( ₁ⁱ , ₂ⁱ), =1,2,..., (3) where M is the number of sensors.

2. Solution of direct problem using the boundary element method

In order to solve the problem (1), (2) for arbitrary assumed values of parameters c1, c2 the Kirchhoff transformation is introduced

µ µ

=

∫

^{λ( )d}

) (

0 T

T

U (4)

and then the governing equations (1) take a form







=

∇

⋅

−

= Γ

∈

=

= Γ

∈

=

∇ Ω

∈

) ( ) ( n ) ( :

) ( ) ( ) ( :

0 ) ( :

2 1

2

x q x U x

q x

T U x U x U x

x U x

b b

b (5)

where (c.f. equation (2))

) ( ) 2 (

1 )

(x c₁T² x c₂T x

U = + (6)

The integral equation corresponding to the problem (5) is the following [1, 2]

Γ

= Γ

+

∫ ∫

Γ

∗ Γ

∗(ξ, ) ( )d (ξ, ) ( )d

) ξ ( ) ξ

( U U x q x Q x U x

B (7)

where ξ is the observation point, B(ξ) is the coefficient from the scope (0, 1], U ^*(ξ,x) is the fundamental solution and for 2D domain oriented in Cartesian coordinate system it is a function of the form

x r

U 1

ln π 2 1 ) , ξ

( =

∗ (8)

where r is the distance between the points ξ and x.

The function

) , ξ ( n ) , ξ

( x U x

Q^∗ = − ⋅∇ ^∗ (9)

is calculated in analytic way and then

r2

π ) 2 , ξ

( d

x

Q^∗ = (10)

where

2 2 2 1 1

1 ξ )cosα ( ξ )cosα

( − + −

= x x

d (11)

and cosα1, cosα2 are the directional cosines of the normal outward vector n.

In numerical realization of the BEM the boundary Γ is divided into N constant boundary elements Γj, j =1,...,N and then the approximation of equation (7) takes a form

∑ ∫

∑ ∫

Γ

∗ Γ =

∗

=

Γ

= Γ +

j j

x Q U x

U q U

B

N

j j N

j j i

i j

i 1

j i 1

d ) , ξ ( d

) , ξ (

(12)

or

j N

j ij j

N

j ij i

iU G q H U

B

1 1

∑

∑

= =

=

+ (13)

where

d ) ξ , (

ⁱ _j

ij U x

G

j

Γ

=

∫

Γ

∗ (14)

while

d ) , ξ (

ⁱ _j

ij Q x

H

j

Γ

=

∫

Γ

∗ (15)

and Bi = B(ξⁱ), Ui = U(ξⁱ), Uj =U(x ^j), qj =q(x ^j). It should be pointed out that for ξⁱ ∈ Γj: Bi = 1/2 (Γj is the constant boundary element), while for ξⁱ ∈ Ω: Bi = 1.

For i =1,2,...,N one obtains the system of N equations of type (13) from which the 'missing' boundary values Uj and qj can be determined.

The values of function Ui, i = N+1, N+2,...,N+L at the internal points are calculated using the formula (c.f. equation (13))

j N

j ij j

N

j ij

i H U G q

U

1 1

∑

∑

= =

−

= (16)

The obtained internal values of function U should be re-counted:

0 ) ( 2

1

2 2

1Ti + c Ti−Ui Ti =

c (17)

One of the roots of this equation corresponds to the searched value of temperature Ti.

3. Sensitivity analysis with respect to c₁ and c₂

Taking into account the dependence (6) the integral equation (7) can be written as follows

 Γ



 



 +

= Γ

+



 



 +

∫

∫

Γ

∗

Γ

∗

d ) ( ) 2 (

1 ) , ξ (

d ) ( ) ξ, ( ξ)

( ξ)

2 ( ) 1 (ξ

2 2 1

2 2 1

x T c x T c x Q

x q x U T

c T

c B

(18)

We differentiate the equation (18) with respect to c1 and then

 Γ



 



 + +

= Γ

+



 



 + +

∫

∫

Γ

∗

Γ

∗

d ) ( )

( ) ( ) 2 ( 1 ) , ξ (

d ) ( ) ξ , ξ) (

( ξ)

( ξ) ( ξ)

2 ( ) 1 (ξ

1 2 1 1 2

1 1

2 1 1 2

x Z c x Z x T c x T x Q

x Q x U Z

c Z T c T

B

(19)

where Z1(x)=∂T(x)/∂c1, Q1(x)=∂q(x)/∂c1.

The approximate form of equation (19) is following



 



 +

= +



 



 +

∑ ∑

=

=

1 2

1 1 1

1

2 λ

2 λ 1

2

1 _{j j j}

N

j j i j

N

j j i i

i i

i T Z G Q H T Z

B (20)

or

2 2

1 1 1

1 1

1 2

1 2

λ 1

λ _{ij j i i}

N

j i i i N

j

j j ij N

j

j

ijQ H Z B Z H T BT

G =

∑

− +

∑

−

∑

= = =

(21)

Next, the equation (18) is differentiated with respect to c2

∫

∫

Γ

∗

Γ

∗

Γ +

+

= Γ +

+ +

]d ) ( Z ) ( ) ( ) ( [ ξ, ) (

d ) ( ) ξ, ( ]

) ξ ( Z ) ξ ( ) ξ ( ) ξ ( [ ) ξ (

2 2 2

1

2 2

2 2

1

x c x T x Z x T c x Q

x Q x U c

T Z

T c B

(22)

where Z2(x)=∂T(x)/∂c2 and Q2(x)=∂q(x)/∂c2. The approximate form of equation (22) is following

(

^{λ 2}

)

^{2 N}

(

^{j λj j2}

)

j ij j

N

j ij i

i i

i T Z G Q H T Z

B + +

∑

=

∑

+

=

=1 1

(23)

or

i N

j

i j ij i

i i N

j

j j ij N

j

j

ijQ H Z B Z H T BT

G

∑ ∑

∑

=

=

=

− +

−

=

1 2 1

2 1

2 λ λ (24)

The boundary conditions are also differentiated with respect to ce, e = 1,2, and then

) 0 ) (

( :

0 ) : (

2 1

∂ =

=∂

∈

∈ =

e e

e

c x x q

Q x

x x Z

Γ Γ

(25)

So, in order to determine the sensitivity functions Zj1 and Zj2 at the boundary nodes x ^j, j= 1,2,...,N two systems of equations (21), (24) should be solved. Next, the values of functions Zi1 and Zi2 at the internal points xⁱ are calculated using again the formulas (21), (24) for i = N+1, N+2,...,N+L. It should be pointed out that the additional problems (21, (24) are coupled with the basic problem (13) because the values of temperatures Tj appearing in (21), (24) should be known.

4. Solution of inverse problem

In order to solve the inverse problem the following least squares criterion is applied [3, 4]

( )

∑

=

−

= ^M

i

di

i T

T c

c S

1

2 2

1, )

( (26)

where ^Ti=T

( )

^x1i,x2i is the calculated temperature T at the point xi for arbitrary assumed values of c1 and c2.

The necessary condition of optimum of function (26) leads to the following system of equations











∂ =

− ∂

∂ =

∂

∂ =

− ∂

∂ =

∂

∑

∑

= =

= =

M

i

c c i i d i M

i

c c i i d i

k ki

c T T c T

S

c T T c T

S

1 2

2

1 1

1

0 )

( 2

0 )

( 2

2 2 1

(27)

where c₁^k,c₂^k for k = 0 are the arbitrary assumed values of parameters c1, c2, while

k

k c

c₁, ₂ for k > 0 result from the previous iteration.

Function Ti is expanded into the Taylor series about known values of c₁^k,c₂^k

( )

c c

(

^{k k}

)

k i k c c k i

i

i c c

c c T c c

T T

T k k 2

1 2 2

1 1 1

1 ^{1 1 2 2}

∂ − +∂

∂ − +∂

= ₌ ⁺ ₌ ⁺ (28)

this means

(

^{k k}

)

^ik

(

^{k k}

)

k i k i

i T Z c c Z c c

T = + _{1 1}⁺¹− ₁ + _{2 2}⁺¹− ₂ (29)

where

k

k c c

k i c i

c k i

i c

Z T c

Z T

2 2 1 1

2 2 1

1 ₌ , ₌

∂

= ∂

∂

= ∂ (30)

are the sensitivity coefficients. Putting (29) into (27) one obtains

( ) ( )

[ ]

( ) ( )

[ ]











=

−

− +

− +

=

−

− +

− +

∑

∑

=

+ +

=

+ +

M

j

k i i d k k k i k k k i k i M

i

k i i d k k k i k k k i k i

Z T c c Z c c Z T

Z T c c Z c c Z T

1

2 2

1 2 2 1

1 1 1 1

1 2

1 2 2 1

1 1 1

0

0

(31)

or

( ) ( )

[ ] ( )

( ) ( )

[ ] ( )











−

=

− +

−

−

=

− +

−

∑ ∑

∑ ∑

= =

+ +

= =

+ +

M

j

M

i

k i k i i d k

i k k k i k k k i M

i

M

i

k i k i i d k

i k k k i k k k i

Z T T Z

c c Z c c Z

Z T T Z

c c Z c c Z

1 1

2 2

2 1 2 2 1

1 1 1

1 1

1 1

2 1 2 2 1

1 1 1

(32)

where k=1,...,K is the number of iteration.

This system of equations can be written in the matrix form ) (

) ( )

( )

(Z^{T k}Z^kc^k+1= Z^{T k}Z^kc^k+ Z^{T k} T_d −T^k (33) where





















=

k M k M

k k

k k

k

Z Z

Z Z

Z Z

2 1

22 21

12 11

...

Z ... (34)

and





















=

















=

k M k k

k

M d

d d

d

T T T

T T T

...

, ...

2 1 2

1

T

T (35)

while

,

₁

2 1 1 1

2

1 



 



=



 



= _k^{k k+} _k^k₊⁺

k

c c c

c c

c (36)

The system of equations (33) allows to determine the values of c₁^k+1,c₂^k+1. The iteration process is stopped when the assumed number of iterations is achieved.

5. Example of computations

The square domain of dimensions 0.1x0.1 m has been considered. On the left surface the boundary heat flux qb = −5 ⋅ 10⁶ W/m² has been assumed, on the remaining parts of the boundary the temperature 100^oC has been accepted.

The boundary has been divided into 40 constant boundary elements.

At first the direct problem has been solved under the assumption that 116

. 386 05427

. 0 ) (

λ T =− T+ (37)

In Figure 1 the position of isotherms 120, 140,..., etc. is shown.

Next, the inverse problem has been considered. It is assumed that the temperatures at four internal nodes (0.02, 0.05), (0.04, 0.05), (0.06, 0.05), (0.08, 0.05) (Fig. 1) are given.

Fig. 1. Temperature distribution

Fig. 2. Identification of parameters for c₁⁰ =−0.3,c₂⁰=300

Fig. 3. Identification of parameters for c₁⁰=−0.3,c₂⁰=700

Fig. 4. Identification of parameters for c₁⁰ =0,c₂⁰=600

Fig. 5. Identification of parameters for c₁⁰=−0.1,c₂⁰=600 - iteration process is not convergent

In Figures 2-4 the results of identification of parameters c1/c1d, c2/c2d (c1d = −0.05426, c2d = 386.116 denote the real values of these parameters) for different initial values c₁⁰, c⁰₂ are shown. It is visible that the iteration process is convergent and after 5 iterations the real values of searched parameters are obtained.

It should be pointed out that iteration process is not always convergent, for example for c₁⁰=0.1 and c⁰₂ =600 - c.f. Figure 5. So, the proper choice of initial values of identified parameters allows the convergence of iteration procedure.

Summing up, the algorithm presented allows to identify the unknown parameter as the temperature dependent function, and it is the main advantage of the approach discussed.

Acknowledgement

This paper has been sponsored by KBN (Grant No 3 T08B 004 28).

References

[1] Brebbia C.A., Dominguez J., Boundary elements, an introductory course, Computational Mechanics Publications, McGraw-Hill Book Company, London 1992.

[2] Majchrzak E., Metoda elementów brzegowych w przepływie ciepła, Wyd. Politechniki Często- chowskiej, Częstochowa 2001.

[3] Kurpisz K., Nowak A.J., Inverse thermal problems, Computational Mechanics Publications, Southampton, Boston 1995, 259-298.

[4] Freus K., Metoda elementów brzegowych w zadaniach odwrotnych ustalonego przepływu ciepła, Praca doktorska, Politechnika Częstochowska, Częstochowa 2002.