Application of boundary element method to shape sensitivity analysis

(1)

APPLICATION OF BOUNDARY ELEMENT METHOD TO SHAPE SENSITIVITY ANALYSIS

Ewa Majchrzak ^1,2, Mirosław Dziewoński ¹, Sebastian Freus ²

1 Silesian University of Technology, Gliwice

2 Czestochowa University of Technology, Czestochowa

Abstract. In the paper the boundary element method to shape sensitivity analysis is applied.

The Laplace equation is analyzed and the aim of investigations is to estimate the changes of temperature in the 2D domain due to the change of local geometry of the boundary. Here the implicit differentiation method of shape sensitivity analysis is used. In the final part of the paper the example of numerical computations is shown.

1. Boundary element method for the Laplace equation

We consider the Laplace equation

( )

²

( )

²

( )

2 2

, ,

, : T x y T x y 0

x y x y

∂ ∂

∈ Ω λ + λ =

∂ ∂ (1)

where λ [W/mK] is the thermal conductivity, T is the temperature, x, y are the geometrical co-ordinates. The equation (1) is supplemented by boundary conditions:

( ) ( )

( ) ( ) ( )

1

2

, : ,

, : , ,

b

x y T x y T

x y q x y T x y q

∈ Γ =

∈ Γ = −λn⋅ ∇ =

(2)

where T_b is the known boundary temperature, q_b is the known boundary heat flux.

In the case considered the boundary integral equation is following [l, 2]

( ) ( ) ( ) ( ) ( )

( ) ( )

*

, : , , , , , , d

, , , , d

B T T x x y q x y

q x x y T x y

Γ

ξ η ∈ Γ ξ η ξ η + η Γ =

η Γ

∫

and it can be calculated analytically

( )

*

, , , 2

2 q x y d

ξ η = r

π (7)

where

( )

x

( )

y

d= x− ξ n + y− η n (8)

while nx, ny are the directional cosines of the normal outward vector n.

In numerical realization of the BEM the boundary is divided into N boundary elements and integrals appearing in equation (3) are substituted by the sums of integrals over these elements

( ) ( ) ( ) ( )

( ) ( )

*

1

*

1

, , , , , , d

, , , , d

j

N

i i i i i i j

j N

i i j

j

B T q x y T x y

T x y q x y

= Γ

ξ η ξ η + ξ η Γ =

ξ η Γ

∑ ∫

(9)

Each point on the linear boundary element Γ_j is expressed as follows (Fig. 1)

(

^,

)

^:

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

 = +

∈ Γ 

= +

 (10)

where

(

xp^j^,yp^j

)

and

(

xk^j^,yk^j

)

correspond to the beginning and end of Γj while

1 1

2 , 2

p k

N − θ N + θ

= = (11)

are the shape functions

(

^{θ ∈ −}

[

^1,1

] )

^.

(3)

Fig. 1. Linear boundary element

It is easy to check that

2 2

d d

d d d d

d d 2 2 2

j j j j

k p k p j

j

x x y y l

x y  −   − 

   

Γ =  θ + θ θ =   +  θ = θ (12)

and

,

j j j j j j

k p y p k

j j x

x y

j j j j

y y l x x l

n n

l l l l

− −

= = = = − (13)

where l_j is the length of the element Γj.

For the linear boundary element Γj, we assume that

( ) ( )

, :

( )

j j

p p k k

j j j

p p k k

T N T N T

x y

q N q N q

 θ = +

∈ Γ 

θ = +

 (14)

The integrals appearing in equation (9) can be written in the form

( ) ( )

* , , , , d ˆ ˆ

j

p j k j

i i j ij p ij k

q x y T x y H T H T

Γ

ξ η Γ = +

∫

⁽¹⁵⁾

and

( ) ( )

* , , , , d

j

p j k j

i i j ij p ij k

T x y q x y G q G q

Γ

ξ η Γ = +

∫

⁽¹⁶⁾

where (c.f. equations (7), (8), (13)):

1

2 1

ˆ 1 d

4

j j j j

x y y x

p

ij p

ij

r l r l

H N

r

−

= θ

π

∫

⁽¹⁷⁾

(4)

1

ln d 4

k j

ij k

ij

G N

r

−

= θ

πλ

∫

⁽²⁰⁾

In equations (17)-(20)

(

^j ^j

) (

² ^j ^j

)

²

( ) ( )

^j ² ^j ²

j k p k p x y

l = x −x + y −y = l + l (21)

and

(

^j ^j

) (

² ^j ^j

)

²

( ) ( )

^j ² ^j ²

ij p p k k i p p k k i x y

r = N x +N x − ξ + N y +N y − η = r + r (22) As is well known, in the final system of algebraic equations the values of temperatures or heat fluxes are connected with the boundary nodes. If the following numeration of the boundary nodes r = l, 2,..., R is accepted then for i = l, 2,.... R one obtains the system of equations (c.f. equation (9))

1 1

R R ˆ

i i ir r ir r

r r

B T G q H T

= =

+

∑

=

∑

⁽²³⁾

where for the single node r being the end of the boundary element Γj and being the beginning of the boundary element Γj+1 (Fig. 2) we have

1

ˆ ˆ ˆ 1

k p

ir ij i j

k p

ir ij i j

G G G

H H H

+

= +

(24)

while for double node r, r +1:

1 1

,

ˆ ˆ , ˆ ˆ

k p

ir ij i r i j

k p

ir ij i r i j

G G G G

H H H H

+ +

= =

= = (25)

The system of equations (23) can be written in the form

1 1

, 1, 2,...,

R R

ir r ir r

r r

G q H T i R

= =

∑ ∑

⁽²⁶⁾

(5)

or

Gq=HT (27)

where

ˆ , ˆ

ir ir

ii i

H i r

H

H B i r

 ≠

= 

− =

 (28)

Fig. 2. Single and double nodes

It should be pointed out that it is convenient to calculate the values H_ii using the formula

1

, 1, 2,...,

R

ii ir

r r i

H H i R

=≠

= −

∑

= ⁽²⁹⁾

Taking into account the boundary conditions (2) the system of equations (27) should be rebuilt to the form A Y = F. The solution of this system allows to deter- mine the “missing” boundary temperatures and heat fluxes. Next, the temperatures at optional set of internal nodes can be calculated using the formula

1 1 1

R R R

i ir r ir r ir

r r r

T H T G q Z

= = =

1

2

, : D 0

D

, : D 0

D x y T

b x y q

b

∈ Γ =

(33)

So, this approach of shape sensitivity analysis is connected with the differentiation of elements of matrices G and H.

Let r

b x

= where

(

x yr^, r

)

is the single boundary node (c.f. Fig. 2). Because

j j1

r k p

x =x =x ⁺ (34)

so after differentiation of G we obtain the following non-zero elements

1 1

, , ,

p k

i j i j

ij ij

j j j j

k k p p

G G

x x x x

+ +

∂ ∂

∂ ∂ ∂ ∂ (35)

and

,

p k

rs rs

r r

G G

∂ ∂

∂ξ ∂ξ (36)

where s = 1, 2,..., j − 1, j + 2,..., N.

Taking into account the formulas (19), (20) one has:

( )

1 1

2

1 1

ln d d

4 4

p j j

ij x j k ir x

p p

j

k j ij ij

G l l N r

N N

x l r r

− −

∂ − δ

= θ − θ

∂ πλ

∫

πλ

∫

⁽³⁷⁾

( )

1 1

2

1 1

ln d d

4 4

k j j

ij x j k ir x

k k

j

k j ij ij

G l l N r

N N

x l r r

− −

∂ − δ

= θ − θ

∂ πλ

∫

πλ

∫

⁽³⁸⁾

(7)

( )

¹

1 1

1

1 1

1 2

1 1 1 1 1

1 1

ln d d

4 4

p j j

p ir x

i j x j

p p

j

k j i j i j

N r

G l l

N N

x l r r

+ +

+

+ − + − +

∂ − δ

= − θ − θ

∂ πλ

∫

πλ

∫

⁽³⁹⁾

( )

¹

1 1

1

1 1

1 2

1 1 1 1 1

1 1

ln d d

4 4

k j j

p ir x

i j x j

k k

j

p j i j i j

N r

G l l

N N

x l r r

+ +

+

+ − + − +

∂ − δ

= − θ − θ

∂ πλ

∫

πλ

∫

⁽⁴⁰⁾

and for s≠ j and s≠ + : j 1

1

2 1

d 4

p s

rs s x

p

r rs

G l r

N r

−

∂ = θ

∂ξ πλ

∫

⁽⁴¹⁾

1

2 1

d 4

k s

rs s x

k

r rs

G l r

N r

−

∂ = θ

∂ξ πλ

∫

⁽⁴²⁾

In similar way, the differentiation of H leads to the following non-zero elements:

( )

²

( )

1

4 1

ˆ 1 2

d 4

j j j j

p

k ir y y ij x ir k

ij j p

k ij

N l r r r N c

H

N

x r

−

 − δ −  + δ −

∂ =   θ

∂ π

∫

⁽⁴³⁾

( )

²

( )

1

4 1

ˆ 1 2

d 4

j j j j

k

k ir y y ij x ir k

ij j k

k ij

N l r r r N c

H

N

x r

−

 − δ −  + δ −

∂ =   θ

∂ π

∫

⁽⁴⁴⁾

( )

¹ ¹ ² ¹

( )

¹

1 1

1

1 4

1 1

ˆ 1 2

d 4

j j j j

p

p ir y y i j x ir p

i j j p

p i j

N l r r r N c

H

N

x r

+ + + +

+ + +

− +

 − δ −  + δ −

∂ =   θ

∂ π

∫

⁽⁴⁵⁾

( )

¹ ¹ ² ¹

( )

¹

1 1

1

1 4

1 1

ˆ 1 2

d 4

j j j j

k

p ir y y i j x ir p

i j j k

p i j

N l r r r N c

H

N

x r

+ + + +

+ + +

− +

 − δ −  + δ −

∂ =   θ

∂ π

∫

⁽⁴⁶⁾

and for s≠ j and s≠ + : j 1

1 2

4 1

ˆ 1 2

d 4

s s s

p

y rs x

rs

p

r rs

l r r c

H N

r

−

− +

∂ = θ

∂ξ π

∫

⁽⁴⁷⁾

1 2

4 1

ˆ 1 2

d 4

s s s

k

y rs x

rs

k

r rs

l r r c

H N

r

−

− +

∂ = θ

∂ξ π

∫

⁽⁴⁸⁾

where

j j j j j

x y y x

c =r l −r l (49)

(8)

tivity equals λ = 1 W/(mK). On the bottom boundary the Neumann condition q_b= −10⁴ W/m² has been assumed, on the remaining parts of the boundary the Dirichlet condition T_b = 500°C has been accepted. The boundary has been divided into 8 linear boundary elements (Fig. 3) and 10 boundary nodes have been distin- guished (two double boundary nodes).

Fig. 3. Discretization

The solution of basic problem (equations (1), (2)) is following

500 10000

676.40 10000

500 10000

500 19122.08

500 324.51

500 , 251.10

500 938.72

500 251.10

500 324.51

500 19122.08

   − 

   − 

   

   − 

   

  − 

=  = 

   

  − 

   

   

T q (50)

(9)

For the shape parameter b = x₆ we have

6

0 0 0 0.009 0.204 0 0 0

0 0 0 0.008 0.225 0 0 0

0 0 0 0 0.240 0 0 0

0 0 0 0 0.299 0 0 0

0.019 0.012 0 0 0.331 0.061 0.038 0.028

0 0 0 0.021 0.254 0 0 0

0 0 0 0.012 0.223 0 0 0

0 0 0 0.012 0.218 0 0 0

0 0 0 0.009 0.204 0 0 0

p

x

 − 

 −



∂ 

= − − − − −

∂ 

 −

 −

 

G



(51)

6

0 0 0 0.015 0.217 0 0 0

0 0 0 0.013 0.234 0 0 0

0 0 0 0 0.237 0 0 0

0 0 0 0 0.288 0 0 0

0.018 0.006 0 0 0.254 0.049 0.036 0.024

0 0 0 0.049 0.331 0 0 0

0 0 0 0.026 0.254 0 0 0

0 0 0 0.023 0.241 0 0 0

0 0 0 0.015 0.217 0 0 0

k

x

 − 

 −



∂ 

= − − − − −

∂ 

 −

 −

 

G



(52)

6

0 0 0 0.773 0.648 0 0 0

0 0 0 1.401 1.126 0 0 0

0 0 0 1.953 1.652 0 0 0

0 0 0 0 3.573 0 0 0

0.457 0.405 1.230 0 0 0 0.710 0.223

0 0 0 0.840 0 0 0 0

0 0 0 0.247 0 0 0 0

0 0 0 0.463 0.405 0 0 0

0 0 0 0.773 0.648 0 0 0

p

x

 − 

 − 

 

 − 

 

 −

∂ 

= − − − − −

∂ 



 −

 

H



(53)

(10)

6 0.498 0.231 1.953 0 0 0 0.563 0.095

0 0 0 0.840 0 0 0 0

0 0 0 0.247 0 0 0 0

0 0 0 0.810 0.868 0 0 0

0 0 0 0.818 0.944 0 0 0

x − − − − −

∂ 

 −

 −

 



 Finally, one obtains the following solution of additional problem (32), (33)

0 0

0.43346 0

0 0

0 1490.11

0 5676.25

D D

0 , 30567.38

D D

0 2117.76

0 384.28

0 140.50

0 15.85

b b

   

   

   − 

   

=  = 

  − 

   

   − 

   

   

T q

(55)

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.