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

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

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 tem- perature 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}

^{,}

### )

^{∈ Ω}

^{:}

^{∇}

^{2}

^{T x y}

### (

^{,}

### )

^{=}

^{0}

^{(1) }

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)

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

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

1

l n 1 d θ 4 π λ

p j

i j p

i j

l

G N

r

−

=

### ∫

^{(5) }

1

1

l n 1 d θ 4 π λ

k j

i j k

i j

l

G N

r

−

=

### ∫

^{(6) }

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

−

−

=

### ∫

^{(7) }

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

−

−

=

### ∫

^{(8) }

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}

^{2}

## ( )

^{j}

^{2}

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)

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}

^{2}

## ( )

^{j}

^{2}

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) (xj^{ }p

, yj^{ }p

), (xjk

, yj^{ }k

) 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-

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) = (xp^{ }j, yp j) then δ_{1} = 1, δ_{2} = 0, δ_{3} = 0, δ_{4} = 0,
- if the observation point (ξi, ηi) = (xk^{ }j

, yk j

) then δ_{1} = 0, δ_{2} = 1, δ_{3} = 0, δ_{4} = 0,

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

, yk j

) then δ_{1} = 0, δ_{2} = 0, δ_{3} = 1, δ_{4} = 0,

- if the observation point (ξ_{i}, η_{i}) not belongs to the circle then δ_{1} = 0, δ_{2} = 0, δ_{3} = 0,
δ_{4} = 1.

We calculate:

## ( ) ^{(} ^{)}

## ( ) ^{(} ^{)}

1 2 3

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

−

∂ ∂ −

=

∂

### ∫

∂ ^{(24) }

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

−

∂ ∂ −

=

∂

### ∫

∂ ^{(25) }

where

## ( )

2

2

4

2

j j j j

x y y x

i j

j 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

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.

It is assumed that λ = 1 W/(mK). On the bottom external boundary the Neumann
condition qb = −10^{4} W/m^{2} 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)

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.