Sensitivity analysis of temperature field with respect to the radiusof internal hole

(1)

SENSITIVITY ANALYSIS OF TEMPERATURE FIELD WITH RESPECT TO THE RADIUS OF INTERNAL HOLE

Ewa Majchrzak ^{1, 2}, Grażyna Kałuża ¹, Marek Paruch ¹

1 Silesian University of Technology, Gliwice

2 Czestochowa University of Technology, Czestochowa

Abstract. The domain (2D problem) with internal hole of radius R is considered.

The temperature field in this domain is described by the Laplace equation supplemented by adequate boundary conditions. The aim of investigations is to estimate the changes of temperature due to change of radius R. In order to solve the problem, the boundary element method is used and the implicit differentiation method of sensitivity analysis is applied.

In the final part of the paper the example of numerical computations is shown.

1. Boundary element method

The steady state temperature field T (x, y) in domain Ω (Fig. 1) is described by the Laplace equation

(

^{x y}^,

)

^{∈ Ω}^: ^∇²^{T x y}

(

^,

)

⁼⁰ ⁽¹⁾

supplemented by adequate boundary conditions.

Fig. 1. Domain considered

Application of the boundary element method [1, 2] leads to the following system of equations

G q = H T (2)

(2)

where vectors T and q contain the values of temperatures T_r and heat fluxes q_r at the boundary nodes (x_r, y_r), r = 1, 2, …, R, G and H are the influence matrices.

If the linear boundary elements are used then for the single node r being the end of the boundary element Γ_j and being the beginning of the boundary element Γ_j+1 we have:

1

ˆ ˆ ˆ 1

k p

i r i j i j

k p

i r i j i j

G G G

H H H

+

= +

(3)

while for double node r, r + 1:

1 1

,

ˆ ˆ , ˆ ˆ

k p

i r i j i r i j

k p

i r i j i r i j

G G G G

H H H H

+ +

= =

(4)

The elements of matrices are calculated using the formulas:

1

l n 1 d θ 4 π λ

p j

i j p

i j

l

G N

r

−

=

∫

⁽⁵⁾

1

l n 1 d θ 4 π λ

k j

i j k

i j

l

G N

r

−

=

∫

⁽⁶⁾

and:

1

2 1

ˆ 1 d θ

4 π

j j j j

x y y x

p

i j p

i j

r l r l

H N

r

−

=

∫

⁽⁷⁾

1

2 1

ˆ 1 d θ

4 π

j j j j

x y y x

k

i j k

i j

r l r l

H N

r

−

=

∫

⁽⁸⁾

where:

1 θ 1 θ

2 , 2

p k

N − N +

= = (9)

are the shape functions (θ ∈[−1, 1]), l_j is the length of the boundary element Γ_j

( )

^j ²

( )

^j ²

j x y

l = l + l (10)

where:

j j j

x k p

j j j

y k p

l x x

l y y

= −

(11)

(3)

In equations (5), (6), (7), (8) rij is the distance between the observation point (ξ_i, η_i) and the point (x_j, y_j) on the linear boundary element Γ_j

( )

^j ²

( )

^j ²

i j x y

r = r + r (12)

while:

ξ η

j j j

x p p k k i

j j j

y p p k k i

r N x N x

r N y N y

= + −

(13)

Fig. 2. Linear boundary element

In formulas (11), (13) (xjp

, yjp

), (xjk

, yjk

) correspond to the beginning and the end of boundary element Γ_j.

2. Shape sensitivity analysis

We assume that radius R of internal hole is the shape parameter (Fig. 1).

The implicit differentiation method [3] of sensitivity analysis starts with the alge- braic system of equations (2). The differentiation of (2) with respect to R leads to the following system of equations

D D D D

DGR + DRq = DHR + DRT

q G T H (14)

This approach of shape sensitivity analysis requires the differentiation of elements of matrices G and H, this means the differentiation of G^p, G^k, H^p, H^k (c.f. equa-

(4)

tions (5), (6), (7), (8)) with respect to the R. Non-zero elements of these matrices are connected with:

1) the integrals over the boundary elements approximating the internal circle, 2) the integrals over the boundary elements approximating the external boundary

of the domain considered for which the observation point (ξi, ηi) belongs to the circle.

In the first case each point on the linear boundary element Γj can be expressed as follows

( ) ( ) ( )

( ) ( )

co s φ co s φ

, :

s i n φ s i n φ

j j

p s p k s k

j j j

p s p k s k

x N x R N x R

x y

y N y R N y R

 = + + +

∈ Γ 

= + + +

 (15)

It is easy to check that (c.f. equation (11)):

( )

co s φ co s φ si n φ s i n φ

j j j j j

x k p k p

j j j j j

y k p k p

l x x R

l y y R

= − = −

(16)

and then

1 ^j ^j

j j x j y

x y

j

l l l

l l

R l R R

 

∂ ∂ ∂

=  + 

∂  ∂ ∂  (17)

where:

c o s φ co s φ

s i n φ s i n φ

j

j j

x

k p

j

y j j

k p

l R l R

∂ = −

∂

= −

∂

(18)

The formulas (13) take the following form:

( ) ( ) ⁽ ⁾ ( )

( )

( ) ( ) ⁽ ⁾ ( )

( )

1 2

3 4

1 2

3 4

δ co s φ δ c o s φ

δ co s φ δ ξ

δ s i n φ δ s i n φ

δ s i n φ δ η

j j j

x p s p k s k

s i i

j j j

y p s p k s k

s i i

r N x R N x R

x R

r N y R N y R

y R

= − + + − + −

+ −

= − + + − + −

+ −

(19)

where:

- if the observation point (ξi, ηi) = (xpj, yp j) then δ₁ = 1, δ₂ = 0, δ₃ = 0, δ₄ = 0, - if the observation point (ξi, ηi) = (xkj

, yk j

) then δ₁ = 0, δ₂ = 1, δ₃ = 0, δ₄ = 0,

(5)

- if the observation point (ξi, ηi) belongs to the circle but (ξi, ηi) ≠ (xpj, ypj) and (ξi, ηi) ≠ (xk j

, yk j

) then δ₁ = 0, δ₂ = 0, δ₃ = 1, δ₄ = 0,

- if the observation point (ξ_i, η_i) not belongs to the circle then δ₁ = 0, δ₂ = 0, δ₃ = 0, δ₄ = 1.

We calculate:

( ) ⁽ ⁾

1 2 3

δ c o s φ δ co s φ δ co s φ

δ s i n φ δ s i n φ δ s i n φ

j

j j

x

p p k k i

j

y j j

p p k k i

r

N N

R r

N N

R

∂ = − + − −

∂

= − + − −

∂

(20)

and then (c.f. equations (5), (6))

1

1 1

1

1 1

l n d θ + 4 π λ

l n 1 d θ 4 π λ

p

i j j

p i j

j

p

i j

G l

N

R R r

l

N

R r

−

∂ ∂

∂ = ∂

 

∂  

∂  

∫

(21)

1

1 1

1

1 1

l n d θ + 4 π λ

l n 1 d θ 4 π λ

k

i j j

k i j

j

k

i j

G l

N

R R r

l

N

R r

−

∂ ∂

∂ = ∂

 

∂  

∂  

∫

(22)

where

2

1 1

ln

j j

y

j x j

x y

i j i j

r r

R r r R R

   ∂ ∂ 

∂   = −  + 

∂    ∂ ∂ 

(23)

In similar way the formulas (7), (8) are differentiated:

1

2 1

ˆ 1

d θ 4π

p j j j j

i j x y y x

p

i j

H r l r l

N

R R r

−

 

∂ ∂ −

=  

∂

∫

∂   ⁽²⁴⁾

1

2 1

ˆ 1

d θ 4π

k j j j j

i j x y y x

k

i j

H r l r l

N

R R r

−

 

∂ ∂ −

=  

∂

∫

∂   ⁽²⁵⁾

(6)

where

( )

2

4

2

j j j j

x y y x

i j

j j

y y

j j j j

x x

y x x y

i j j j

j x j y j j j j

x y x y y x

i j

r l r l

R r

l r

r l

l r l r

R R R R

r r r

r r r l r l

R R

r

 − 

∂   =

∂  

∂ ∂ ∂ ∂ 

+ − −

 

 ∂ ∂ ∂ ∂ 

  −

 ∂ ∂ 

+ −

 

 ∂ ∂ 

 

(26)

In the second case for the integrals over the boundary elements approximating the external boundary of the domain considered for which the observation point (ξ_i, η_i) belongs to the circle we have

(

^,

)

^:

j j

p p k k

j j j

p p k k

x N x N x

x y y N y N y

 = +

∈ Γ  = + (27)

and then

0, 0, 0

j j j

x y

l l l

R R R

∂ ∂ ∂

= = =

∂ ∂ ∂ (28)

while:

( )

i

co s φ si n φ

j j j

x p p k k s

j j j

y p p k k s

r N x N x x R

r N y N y y R

= + − +

= + − + (29)

from which:

co s φ

s i n φ

j x

i j

y

i

r R r

R

∂ = −

∂

= −

∂

(30)

The next part of the algorithm is similar as for the first case.

3. Example of computations

The square of dimensions 0.05×0.05 m is considered - Figure 3. The center of the circle: xs = 0.025, ys = 0.025, radius: R = 0.015 m.

(7)

It is assumed that λ = 1 W/(mK). On the bottom external boundary the Neumann condition qb = −10⁴ W/m² is accepted, on the remaining part of the external boundary the Dirichlet condition T_b1 = 500°C is assumed. Along the circle the constant temperature T_b2 = 700°C is given.

Fig. 3. Discretization

The solution of the basic problem is following

5 0 0 1 0 0 0 0

7 8 2 . 8 3 9 3 1 0 0 0 0

5 0 0 1 0 0 0 0

5 0 0 2 4 9 2 1 . 7 1

5 0 0 7 5 5 . 3 5

5 0 0 8 6 6 . 3 7

5 0 0 1 8 7 2 8 . 0 8

5 0 0 , 8 6 6 . 3 7

5 0 0 1 5 7 5 5 . 3 5

5 0 0 2 4 9 2 1 . 7 1

7 0 0 2 6 3 0 4 . 0 1

7 0 0 7 5 6 7 . 9 5

7 0 0 2 6 3 0 4 . 0 1

7 0 0 2 2 2 0 0 . 6 0

  −

  −

 

  −

 

= =

 

  −

 

  −

 

 

T q

 

 

 

(31)

(8)

while the solution of the additional problem connected with the sensitivity analysis

0 0

6 2 5 . 6 5 1 8 0

0 0

0 1 0 9 0 1 8 . 3 4

0 1 9 4 2 0 9 6 . 9 5

0 1 2 1 6 3 . 0 3

0 2 1 3 6 0 3 8 . 2 9

D D

0 , 1 2 1 6 3 . 0 3

D D

0 1 9 4 2 0 9 6 . 9 5

0 1 0 9 0 1 8 . 3 4

0 1 5 5 5 7 0 3 . 6 5

0 1 4 9 2 9 0 0 . 5 1

0 1 5 5 5 7 0 3 . 6 5

0 8 3 0 2 5 8 . 6 9

R R

  

  

   −

  

  −

 

= =

  −

 

  −

 

  −

 

  −

 

  

T q





 

 

(32)

Using the expansion of function T into the Taylor series

( ) ( )

D

( )

D

T R R T R T R R

+ ∆ = + R ∆ (33)

we can estimate the change of temperature due to the change of R. For example, if

∆R = 0.1R then the temperature at the node 2 (Fig. 3) changes from 782.8393 to 782.8393 + 625.6518 ⋅ 0.0015 = 783.7778.

References

[1] Brebbia C.A., Domingues J., Boundary elements, an introductory course, CMP, McGraw-Hill Book Company, London 1992.

[2] Majchrzak E., Boundary element method in heat transfer, Publ. of the Techn. Univ. of Czest., Czestochowa 2001 (in Polish).

[3] Burczyński T., Sensitivity analysis, optimization and inverse problems, (in:) Boundary element advances in solid mechanics, Springer-Verlag, Wine, New York 2004.