Design sensitivity analysis and optimization of interface shape for zoned-inhomogeneous thermal conduction problems using boundary integral formulation

ARTICLE IN PRESS Engineering Analysis with Boundary Elements 34 (2010) 825–833

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Design sensitivity analysis and optimization of interface shape for zoned-inhomogeneous thermal conduction problems using boundary integral formulation Boo Youn Lee n Department of Mechanical & Automotive Engineering, Keimyung University, 1000 Shindang-dong, Dalseo-gu, Daegu 704-701, Republic of Korea

a r t i c l e in fo

abstract

Article history: Received 18 January 2010 Accepted 5 May 2010 Available online 3 June 2010

A generalized formulation of the shape design sensitivity analysis for two-dimensional steady-state thermal conduction problem as applied to zoned-inhomogeneous solids is presented using the boundary integral and the adjoint variable method. Shape variation of the external and zone-interface boundary is considered. Through an analytical example, it is proved that the derived sensitivity formula coincides with the analytic solution. In numerical implementation, the primal and adjoint problems are solved by the boundary element method. Shape sensitivity is numerically analyzed for a compound cylinder, a thermal diffuser and a cooling fin problem, and its accuracy is compared with that by numerical differentiation. The sensitivity formula derived is incorporated to a nonlinear programming algorithm and optimum shapes are found for the thermal diffuser and the cooling fin problem. & 2010 Elsevier Ltd. All rights reserved.

Keywords: Shape design sensitivity analysis(SDSA) Shape optimization Thermal conduction Zoned-inhomogeneous solids Boundary element method

1. Introduction The shape optimization problem is to find a continuum shape minimizing an objective function under prescribed constraints. The shape design sensitivity analysis (SDSA) is required to predict the sensitivity of state variables or performance functionals in numerical implementation of shape optimization algorithms. Accurate prediction of the shape sensitivity is very important because inaccurate sensitivity may cause much computational time and moreover lead to undesirable solution in numerical optimization. Hence, the SDSA has been one of intensive subjects to researchers of the shape optimization. Nowadays, the boundary element method (BEM) or the boundary integral equation (BIE) has especially received considerable attention in the area of the SDSA. One can refer Burczynski [1], Lee [2], Choi et al. [3] and Li et al. [4] for a survey of research on the BEM-based SDSA. The methods can be classified into a semi-analytical and an analytical approach. In the semianalytical approach, the state equations are first discretized into boundary elements and next the derivatives of the system matrices are derived analytically. In the analytical approach, the material derivatives [5] are taken to the state equations and sensitivity formulas are analytically derived on the continuum basis, next discretization followed only for the numerical implementation. The continuum approach can again be classified

n

Tel.: + 82 53 580 5922; fax: + 82 53 580 5115. E-mail address: [email protected]

0955-7997/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2010.05.009

into the adjoint variable method (AVM) and the direct differentiation method (DDM). In the AVM, the shape sensitivity formulas are derived by introducing adjoint systems that are determined depending on the performance functional, whereas, in the DDM, the sensitivity of the state variables is directly calculated by solving the BIE for the design derivatives. While most of the research on the SDSA has dealt with the elasticity problem, papers focusing on the SDSA of solids in thermal environment are not so many. Most of the published papers for the SDSA of thermal conduction problems belong to the analytical approach. Meric [6–8] has derived shape sensitivity formulas for thermoelastic and thermal conduction problems using the Lagrangian multiplier technique, the material derivative concept and the AVM. Park and Yoo [9] have applied the method of variational formulation [5] to the thermal conduction problem and derived a sensitivity formula, transforming the variational equation to an equivalent BIE. The SDSA based on the BIE formulation has been presented for thermoelasticity problem (Lee and Kwak [10,11]) and for thermal conduction problem (Lee et al. [12,13] and Lee [14]), using the AVM and the DDM. Aithal and Saigal [15] have presented the SDSA of Dirac delta type performance functional for the thermal conduction problem using the AVM. Sluzalec and Kleiber [16] have presented a shape sensitivity formula for a nonlinear thermal conduction problem with temperature-dependent thermal conductivity using the AVM. Kane and Wang [17,18] presented the semi-analytical approach for nonlinear thermal conduction problems with nonlinear boundary conditions and temperature-dependent thermal conductivity.

ARTICLE IN PRESS 826

B.Y. Lee / Engineering Analysis with Boundary Elements 34 (2010) 825–833

Most of the previous BEM-based researches for the SDSA of thermal conducting solids are concerned with homogeneous problems. Present paper deals with the SDSA of inhomogeneous thermal conduction problems with zoned-interfaces of different materials. A cooling fin, a molding die, a composite tube and so on, consisting of several materials with different thermal conductivities can be examples of such problems. Choi et al. [3] presented the SDSA of an interface problem as applied to implant design in dentistry. Their method uses the BIE formulation and is confined to the elasticity problem. Dems and Mroz [19] have presented the SDSA of the two-dimensional thermal conduction problems with zoned-interfaces, and they derived a general shape sensitivity formula for varying interfaces using the variational formulation based on the finite element method. Recently, Gao and He [20] presented a finite-difference approach called as the complex-variable-differentiation method for the shape sensitivity of multi-region heat conduction using the BEM. The method of the SDSA in present paper is basically different from Dems and Mroz [19] and Gao and He [20] in that it is an analytic approach based on the BIE formulation. The material derivative concept is used to represent the shape variation of the external and interface boundaries. The shape sensitivity formula is derived using the boundary integral identity [12,13] and employing the AVM. The theoretical formulation is validated with an analytical example of a composite wall problem. Shape sensitivity is numerically analyzed and its accuracy is verified for three numerical examples: a compound cylinder, a thermal diffuser and a cooling fin problem. As application to numerical optimization, the sensitivity formula is incorporated to a nonlinear programming algorithm and optimum shapes are found for the thermal diffuser problem and the cooling fin problem.

flow rate per unit area is expressed as qi ¼ kT,i

ð1Þ

where k is the thermal conductivity of the material. The indices in the tensor notations have a range of 1 and 2. For a case without heat source, the steady-state thermal equilibrium equation is given by qi,i ¼ 0 in O

ð2Þ

A general thermal boundary condition is applied on the external boundary G ¼ GT[Gq[Gh as T ¼ T0 at GT q  qi ni ¼ q0 at Gq q  qi ni ¼ hðTT1 Þ at Gh

ð3Þ

where q is the normal heat flux, ni the outward unit normal vector on the boundary. Temperature T0 is prescribed at GT, normal heat flux q0 at Gq and constant heat transfer coefficient h and surrounding temperature TN at Gh. At the interface GI, continuity of the temperature and the normal heat flux should be imposed as T ð1Þ ¼ T ð2Þ ¼ T c qð1Þ ¼ qð2Þ ¼ qc

at GI

ð4Þ

Introducing an arbitrary system of temperature Tn and normal heat flux qn and using the procedure in Lee et al. [12,13], one can derive the boundary integral identity (BII) applicable to the zoned-inhomogeneous thermal conduction problem as Z ðT ð1Þ qð1Þ qð1Þ T ð1Þ Þds ¼ 0 Gð1Þ [GI E

Z

ðT ð2Þ qð2Þ qð2Þ T ð2Þ Þds ¼ 0

ð5Þ

Gð2Þ [GI E

2. Method of the SDSA Consider a two-dimensional isotropic and inhomogeneous thermal conducting solid body O ¼ O(1)[O(2) with zoned-interface boundary GI, as shown in Fig. 1. The superscript (1) or (2) denotes a zone number identifying a material type. External boundary of ð2Þ ð1Þ ð2Þ O is represented as G ¼ Gð1Þ E [ GE , where GE and GE represents external boundary of each zone. The position in the solid body will be denoted by x when necessary. If the temperature is denoted by T(x), then the heat flux vector representing the heat

(2)

ΓE (1 )

ΓE

ΓI Ω

Ω

(2)

(1 )

y Ω =Ω

(1 )

(1 )

∪ Ω

(2)

(2)

Γ = ΓE ∪ Γ E

x Fig. 1. A two-dimensional thermal conducting solid body with zoned-interface.

where ds represents integration with respect to x on the boundary. Above equations correspond to Green’s second identity that holds between two thermal equilibrium states: one with T(i) and q(i), and the other with T(i)n and q(i)n. Now, consider a general performance functional arising in the shape optimization in the following form Z F ¼ cðT,qÞds ð6Þ G

We define the performance functional F on the external boundary G, taking into consideration that most of the shape optimization problems for the thermal conducting solids are concerned with the state variables on the external boundary. Shape variation in the shape optimization problem is defined as a change of the shape of the external boundary G and the interface boundary GI. Using the material derivative concept [5], we can define the shape variation using a velocity vector V  fV1 ,V2 g and a time-like parameter t as shown in Fig. 2. Taking the material derivative to Eq. (6), variation of F can be expressed as Z  Z  : : dF ¼ cT T þ cq q ds þ cVj,s sj ds G G Z Z Z Z : : : cq qdsþ cT T dsþ ðcT þ hcq ÞT ds þ cVj,s sj ds ð7Þ ¼ GT

Gq

Gh

G

_ where cT and cq denote partial derivatives, and T_ and qthe material derivatives, respectively. Vj represents the design velocity vector, sj a unit tangential vector on the boundary, and a subscript (,s) the differentiation in the tangential direction on the boundary. On the other hand, taking the material derivative to each BII in Eq. (5), summing the resulting two equations, and

ARTICLE IN PRESS B.Y. Lee / Engineering Analysis with Boundary Elements 34 (2010) 825–833

Z þ

ΓE( 2)

ΓE(1)

G

Ω (1) x

τV(x) x’

Fig. 2. Shape variation.

ð8Þ

ð13Þ

Taking the material derivative of Fk, we obtain Z  Z xT _ xq x dL dFk ¼ Vj,s sj dsFk k T þ q_ dsþ Lk Lk Gk Lk Gk Lk

ð16Þ

where dLk represents variation of the length of Gk due to shape variation. The boundary condition of the adjoint problem as described in Eq. (11) uses the following substitution:

   RðVÞ ¼ qc sj Vj T,sð1Þ þ T,sð2Þ Vj,s sj ðT ð1Þ T ð2Þ Þ þ nj Vj ðqð1Þ =kð1Þ      þqð2Þ =kð2Þ Þ T,sc nj Vj kð1Þ T,sð1Þ þ kð2Þ T,sð2Þ þsj Vj ðqð1Þ þ qð2Þ Þ ð9Þ and the thermal conductivity of each zone is defined as k(1) and k(2). Then substituting Eq. (8) into Eq. (7) gives Z Z Z : : : dF ¼ ðcq T  Þqds þ ðcT þq ÞT dsþ fcT þ hðcq T  Þ þq gT ds GT

Gq

Gh

Z n: o :  qc ðT ð1Þ T ð2Þ ÞT c ðqð1Þ þ qð2Þ Þ ds ZGI þ fðqsj þ kT,s nj ÞVj T,s þ ðqnj =kT,s sj ÞVj q qT  Vj,s sj gds Z ZG þ cVj,s sj ds þ RðVÞds

at GT

q ¼ hfT  ðcq þ cT =hÞg

at Gh

at Gq

at GI

xT

,

Lk

cq ¼

Z

þ

Gk

ð11Þ

ð12Þ

Substituting Eq. (11) and (12) into Eq. (10), we finally obtain a sensitivity formula as follows: Z dF ¼ fðqsj þkT,s nj ÞVj T,s þðqnj =kT,s sj ÞVj q qT  Vj,s sj gds

ð17Þ

Lk

G

þ

ð10Þ

xq

Then the sensitivity formula of Fk is finally expressed as follows: Z dFk ¼ fðqsj þ kT,s nj ÞVj T,s þ ðqnj =kT,s sj ÞVj q qT  Vj,s sj gds

GI

GI

T  ¼ cq q ¼ cT

cT ¼

Z

Now, in order to remove the material derivatives of the state _ we define an adjoint system with the boundary variables, T_ and q, condition given in Eq. (11) and continuity at the interface given in Eq. (12).

G

 1 1 c c ð1Þ ð2Þ T  T ðk þ k Þ nj Vj ,s ,s kð1Þ kð2Þ

Gk

where

qð1Þ ¼ qð2Þ ¼ qc

GI



where Lk denotes length of the boundary segment Gk given by Z ds ð15Þ Lk ¼

GI

T ð1Þ ¼ T ð2Þ ¼ T c

qc qc

Eq. (13) represents sensitivity of the general performance functional F with respect to shape variation of the external boundary G and the interface GI. Since present formulation is based on the BIE formulation, we can naturally use the BEM to solve the primary problem with the boundary condition and continuity, Eqs. (3) and (4), and the adjoint problem with Eqs. (11) and (12). And then transforming the shape variation into the design velocity vector V, we can calculate the shape sensitivity by numerical integration of Eq. (13). As a specific case of the performance functional, we will derive a sensitivity formula of a functional Fk representing an average value of a function x(T,q) at the boundary segment Gk as follows: R xðT,qÞds Fk ¼ Gk ð14Þ Lk

ΓI

G

Z

2qc T,sc sj Vj ds

Ω ( 2)

applying the continuity at the interface, we obtain Z Z c _  T_ q Þdsþ ðqT fq_ c ðT ð1Þ T ð2Þ ÞT_ ðqð1Þ þ qð2Þ Þgds G GI Z ¼ fðqsj þkT,s nj ÞVj T,s þðqnj =kT,s sj ÞVj q qT  Vj,s sj gds GZ þ RðVÞds

cVj,s sj ds þ

827

x Lk

qc qc



  1 1  ð2Þ T,sc T,sc ðkð1Þ þ kð2Þ Þ nj Vj 2qc T,sc sj Vj ds ð1Þ k k

Vj,s sj dsFk

dLk Lk

ð18Þ

If we need to define the performance functional as average temperature or average heat flux at a boundary segment, the sensitivity formula of Eq. (18) can be efficiently used.

3. Analytical example In this section the sensitivity formulation described in the previous section is validated with an analytical example of a composite wall problem as shown in Fig. 3. Upper and lower boundaries of the wall are insulated. As a boundary condition of the wall, temperature T0 is prescribed at G1, and a convection boundary condition with constant heat transfer coefficient h and surrounding temperature TN at G2. Thermal conductivity of O1 and O2 is defined as k1 and k2, respectively. We want to derive the sensitivity of the temperature and heat flux at G2 in terms of the location of the interface GI, i.e., a. Analytical solution of the problem can be easily obtained from the one-dimensional thermal equilibrium equation using the left and right boundary conditions and continuity at the interface.

ARTICLE IN PRESS 828

B.Y. Lee / Engineering Analysis with Boundary Elements 34 (2010) 825–833

y

Taking into consideration that the design velocity V1 ¼ da at GI and V1 ¼0 elsewhere, and substituting Eq. (24) into Eq. (13), we obtain

  1 1 dFT ¼ qc qc  da ð25Þ k1 k2

Γ1

T=T0

Ω1

Ω2

ΓI

Γ2

h, Τ∞ x

a

Substitution of qc in Eq. (19) and qcn in Eq. (24) into Eq. (25) gives the same result of the sensitivity of the temperature as given in Eq. (21). It is omitted but can be easily proved that sensitivity of the heat flux by the present method of the SDSA gives the same result as the analytical solution in Eq. (21).

l 4. Numerical examples

Fig. 3. An analytical example: compound wall problem.

Solution of temperature and heat flux is expressed as   k0 aT1 þ kh2 þla T0 Tc ¼ C   k2 T þ l þ ðk0 1Þa T1 0 T9G2 ¼ h C k2 ðT0 T1 Þ c q ¼ q9G2 ¼ C

ð19Þ

where k2 k1 k2 þ l þ ðk0 1Þa C¼ h k0 ¼

ð20Þ

   Differentiating the temperature and heat flux at G2, T G2 and qG2 in Eq. (19), with respect to a, we obtain the sensitivity as follows: dðT9G2 Þ da dðq9G2 Þ da

¼

ð1k0 Þk2 ðT0 T1 Þ hC 2

¼

ð1k0 Þk2 ðT0 T1 Þ C2

ð21Þ

Next, using the method of the SDSA presented in the previous section, sensitivity of the following two functionals are derived: Z FT ¼ Tds

Fq ¼

Z

Three example problems with interface shape of zonedinhomogeneous thermal conducting solid are considered to illustrate the numerical implementation. For a compound cylinder, a thermal diffuser, and a cooling fin problem, shape sensitivity is numerically analyzed using the derived formula and its accuracy is compared with that by the method of finite difference. In the method of finite differences, 0.1% change of the design variables is used for the numerical differentiation. A BEM program for two-dimensional zoned-inhomogeneous thermal conduction problems is developed to solve the primal and adjoint systems using quadratic boundary elements. As application to shape optimization, the sensitivity formula is incorporated in the optimization algorithm by Arora [21] and optimum shapes are found for the thermal diffuser problem and the cooling fin problem. 4.1. Compound cylinder problem Consider a compound cylinder with triple layers of which thermal conductivities are k1, k2 and k3, as shown in Fig. 4. It is assumed as a two-dimensional thermal conduction problem. Surfaces of the cylinder are exposed to convection boundary conditions. At inner surface, heat transfer coefficient h1 and surrounding temperature and TN1 are prescribed, and at outer surface, h2 and TN2. Numerical values are set as a ¼0.5, b¼0.6, c ¼0.7, d ¼0.8 m; h1 ¼ 100, h2 ¼20 W/m2 1C; TN1 ¼100, TN2 ¼20 1C; k1 ¼50, k2 ¼30, k3 ¼20 W/m 1C. Considering symmetry, the section of the cylinder is modeled by 70 quadratic boundary elements as shown in Fig. 5. The dots in Fig. 5 represent end-nodes of the quadratic elements.

G2

y qds

ð22Þ

G2

where height of the wall, i.e., length of G2 is assumed as 1. External boundary condition of the adjoint problem for the temperature functional can be expressed, from Eq. (11), as  T  G1 ¼ 0   q G ¼ hðT  G 1=hÞ ð23Þ 2

d

c b x

2

a Solution of the adjoint problem can be obtained from the equilibrium equation using the boundary condition Eq. (23) and continuity Eq. (12), as follows: k0 a T ¼ hC k2 qc ¼  hC

h1,T∞ 1 k1

c

h2,T∞ 2 ð24Þ

k2 k3

Fig. 4. A compound cylinder problem.

ARTICLE IN PRESS B.Y. Lee / Engineering Analysis with Boundary Elements 34 (2010) 825–833

Design boundary is defined as the interface G9, and hence radius b is the design variable. The performance functionals of the temperature and heat flux are defined as follows: R Tds F T ¼ Gk , Gk A G2 Lk R ð26Þ qds F q ¼ Gk , Gk A G4 Lk where Gk coincides to each boundary element on G2 or G4. The temperature functional FT is defined on element number 6–10, and the heat flux functional Fq on element number 16–25. Predicted sensitivities of the temperature and heat flux functionals are listed in Table 1 and are compared with those from the finite differences. The ratio in Table 1 is determined by dividing the sensitivity prediction by the finite difference. The maximum error with respect to the finite difference result is 0.28% for the temperature functional and 0.04% for the heat flux functional. It is shown that very accurate results are obtained. 4.2. Thermal diffuser problem Consider a shape optimization problem for a thermal diffuser with two zones of different materials, as shown in Fig. 6. Thermal

829

conductivities of the two regions, k1 and k2, are 50 and 100 W/ m 1C, respectively. A uniform inward heat flux of q0 ¼ 100 kW/m2 is given at the heat source surface G1 and G2, and a convection boundary condition with heat transfer coefficient h¼500 W/m2 1C and surrounding temperature TN ¼20 1C at the heat discharge surface G4 and G5. The lateral surfaces G6 and G3 are insulated. Design variables are defined as y-coordinates of the two end points of the interface G7, a and b, and design boundary G7 is represented by a line connecting the two end points. The thermal diffuser is modeled by 50 quadratic boundary elements as shown in Fig. 7. Design sensitivities are analyzed for the performance functionals of the temperature and heat flux as follows: R Tds F T ¼ Gk , Gk A G1 [ G2 Lk R ð27Þ qds F q ¼ Gk , Gk A G4 [ G5 Lk where the temperature functional FT is defined on element number 21–30 at the heat source surface G1 and G2, and the heat flux functional Fq on element number 1–10 at the heat discharge surface G4 and G5. For a case of {a, b} ¼{0.025, 0.1} m, results of analyzed sensitivities of the temperature and heat flux functionals are listed in Table 2 and are compared with those from the finite differences. The maximum error with respect to the finite

y

Γ5 Γ6

Γ Γ

k2

Γ7 Γ1 Γ2

Γ4

b

k1

0.05 m a

0.2 m

Γ3

x

Γ

0.2 m Fig. 5. Boundary element model and element numbers of the compound cylinder.

Fig. 6. A thermal diffuser problem.

Table 1 Sensitivity of the compound cylinder problem.

FT

Fq

Element number

Functional values

Finite difference

Sensitivity prediction

Ratio (%)

6 7 8 9 10

0.802404D +02 0.793622D+ 02 0.785114D+ 02 0.776863D+ 02 0.768856D+ 02

 0.159111D  01  0.100826D  01  0.461547D  02 0.522826D  03 0.536166D  02

 0.159023D  01  0.100840D  01  0.461420D  02 0.524311D 03 0.536364D 02

99.94 100.01 99.97 100.28 100.04

16 17 18 19 20 21 22 23 24 25

0.102081D +04 0.102081D +04 0.102081D +04 0.102081D +04 0.102081D +04 0.102081D +04 0.102081D +04 0.102081D +04 0.102081D +04 0.102081D +04

0.138892D+ 00 0.138891D+ 00 0.138891D+ 00 0.138891D+ 00 0.138891D+ 00 0.138891D+ 00 0.138891D+ 00 0.138891D+ 00 0.138891D+ 00 0.138892D+ 00

0.138952D+ 00 0.138929D+ 00 0.138937D+ 00 0.138938D+ 00 0.138939D+ 00 0.138939D+ 00 0.138938D+ 00 0.138937D+ 00 0.138929D+ 00 0.138952D+ 00

100.04 100.03 100.03 100.03 100.03 100.03 100.03 100.03 100.03 100.04

ARTICLE IN PRESS 830

B.Y. Lee / Engineering Analysis with Boundary Elements 34 (2010) 825–833

difference result is 2.15% on the element number 25 for the temperature functional, and 5.97% on the element number 5 for the heat flux functional. The elements 5 and 25, of which sensitivity errors are the maximum, are located at corners where one zone interfaces with the other. Generally applying the BEM to such multiple-zoned thermal problems, supplemental equations are required to deal with multiple heat fluxes defined on the corner nodes. In dealing with zoned-inhomogeneous heat conduction problems, discontinuity of heat flux at corner points shared by different materials is a big issue [22]. In the present research, the linear interpolation method [23] for the temperature field in a triangle, of which vertices are defined by a corner node and two adjacent nodes, is used as the supplemental equation, and hence heat flux in the triangle is approximated to be constant. The error in the elements adjacent to the corners is caused from the BEM implementation, not from the sensitivity formula, which

may be improved in the future. Nevertheless, the sensitivity results listed in Table 2 appear fairly accurate overall. The shape optimization problem is to determine y-coordinates of the two end points of the interface G7, a and b, minimizing the difference of the maximum and the minimum heat fluxes at the heat discharge surface G4 and G5 as follows: Find fa,bg

Fqmax Fqmin

minimizing F ¼ R G qds where Fq ¼ k , Lk

q0

Gk A G4 [ G5

ð28Þ

Upper and lower bounds of the design variables are set as follows: 0:01 ra r0:04 0:01 rb r 0:19

ð29Þ

Initial design variables are set as {a, b} ¼{0.025, 0.1} m and numerical optimization is executed. The optimal design is obtained in 54 iterations from the initial design and optimal design variables are {a, b} ¼{0.0150, 0.1163} m. The initial and

Γ Initial Optimal

Γ

Γ Γ

Fig. 7. Boundary element model and element numbers of the thermal diffuser.

Fig. 8. Initial and optimal shapes of the thermal diffuser.

Table 2 Sensitivity of the thermal diffuser problem.

FT

Fq

Element number

Functional values

Finite difference

Sensitivity prediction

Ratio (%)

21 22 23 24 25 26 27 28 29 30

0.203642D + 03 0.206987D + 03 0.210310D + 03 0.213834D+ 03 0.217994D+ 03 0.224570D + 03 0.229690D + 03 0.232809D + 03 0.234710D + 03 0.235625D+ 03

0.413790D  01 0.395897D 01 0.384736D 01 0.377070D  01 0.367785D 01 0.399082D  01 0.452026D  01 0.484232D 01 0.503951D  01 0.513258D 01

0.415196D  01 0.394159D  01 0.381349D  01 0.370363D  01 0.375701D  01 0.404999D  01 0.451647D  01 0.483904D  01 0.503190D  01 0.515247D  01

100.34 99.56 99.12 98.22 102.15 101.48 99.92 99.93 99.85 100.39

1 2 3 4 5 6 7 8 9 10

0.256992D+ 05 0.256090D + 05 0.254736D+ 05 0.254112D+ 05 0.257984D+ 05 0.265079D + 05 0.258174D+ 05 0.246614D+ 05 0.232777D+ 05 0.217697D+ 05

0.242296D+ 01 0.184169D+ 01 0.748015D + 00  0.710112D + 00  0.218913D+ 01  0.300645D +01  0.281716D+ 01  0.198384D+ 01  0.812469D+ 00 0.551380D + 00

0.242708D + 01 0.184800D + 01 0.757597D +00  0.686964D +00  0.231986D +01  0.298101D + 01  0.279569D +01  0.196957D +01  0.801041D + 00 0.564108D + 00

100.17 100.34 101.28 96.74 105.97 99.15 99.24 99.28 98.59 102.31

ARTICLE IN PRESS B.Y. Lee / Engineering Analysis with Boundary Elements 34 (2010) 825–833

27000 26000 25000 Heat flux

optimal shapes are shown in Fig. 8. Values of the objective function F for the initial and optimal designs are 0.04738 and 0.04178, respectively. Fig. 9(a) and (b) represent temperature distribution for the initial and optimal designs, respectively. Distribution of the heat flux at the heat discharge surface is represented in Fig. 10, showing that the difference of the maximum and the minimum heat fluxes is decreased at the optimal design in comparison with the initial design.

831

24000 23000 C

4.3. Cooling fin problem

22000 Consider a shape optimization problem for a cooling fin with two zones of different materials, as shown in Fig. 11. Thermal conductivities of the wall and fin, k1 and k2, are 50 and 100 W/m 1C, respectively. Convection boundary conditions are prescribed as

B 21000

Initial design Optimal design

A

20000 A

B

C

Position

60

Fig. 10. Distribution of heat flux at the heat discharge surface of the thermal diffuser (unit: W/m2).

80 100 Γ

120 Γ

160 Γ

Γ Γ

200

Fig. 11. A cooling fin problem.

220

180

140

60

Γ

80 Γ

100 Γ

120 Fig. 12. Boundary element model and element numbers of the cooling fin.

160 200

220

180

140

Fig. 9. Temperature distribution of the thermal diffuser at the initial and optimal designs(unit: 1C). (a) Initial design. (b) Optimal design.

h1 ¼200 W/m2 1C, TN1 ¼100 1C at the heat source surface G5, and h2 ¼20 W/m2 oC, TN2 ¼20 oC at the heat discharge surface G1, G2 and G3. If fins are located in the same intervals along the long wall, then hatched upper and lower surfaces shown in Fig. 11 can be assumed insulated considering symmetry. Three design variables are defined as y-coordinates of the two end points of the lateral surface G2, a and b, and length of the fin, c, the lateral surface G2 being represented by a line connecting the two end points. The wall and fin are modeled by 65 quadratic boundary elements as shown in Fig. 12.

ARTICLE IN PRESS 832

B.Y. Lee / Engineering Analysis with Boundary Elements 34 (2010) 825–833

Table 3 Sensitivity of the cooling fin problem. Element number

Functional values

Finite difference

Sensitivity prediction

Ratio (%)

21 23 25 27 29 31 33 35 37 39 41 43 45

0.105003D + 02 0.104910D + 02 0.104684D + 02 0.213587D +02 0.216434D +02 0.220068D + 02 0.224358D +02 0.229225D +02 0.234667D +02 0.240947D + 02 0.248278D +02 0.253594D +02 0.255163D +02

0.824064D  02 0.822404D  02 0.818373D  02 0.168303D  01 0.173375D  01 0.179865D  01 0.187549D  01 0.196293D  01 0.206081D  01 0.217290D  01  0.269021D  01  0.277417D  01  0.279893D  01

0.811289D  02 0.809757D  02 0.805612D  02 0.165725D  01 0.170756D  01 0.177198D  01 0.184822D  01 0.193481D  01 0.203095D  01 0.213566D  01  0.258417D  01  0.274961D  01  0.277943D  01

98.45 98.46 98.44 98.47 98.49 98.52 98.55 98.57 98.55 98.29 96.06 99.11 99.30

Fig. 13. Initial and optimal shapes of the cooling fin.

Design sensitivities are analyzed for the heat flux functional as follows: Z Fq ¼ qds, Gk A G1 [ G2 [ G2 ð30Þ Gk

where the boundary segment Gk is defined as element number 21–45 at the heat discharge surface G1, G2 and G3. For a case of {a, b, c}¼ {0.10, 0.05, 0.30} m, results of analyzed sensitivity of the heat flux functional are listed in Table 3 and are compared with those from the finite differences. The sensitivity results appear fairly accurate overall. The maximum error with respect to the finite difference result is 3.94% on the element number 41 adjacent to the zone corner. The shape optimization problem is to determine the dimension of the fin {a, b, c}, maximizing an integral of the heat flux along the heat discharge surface G1, G2 and G3 as follows: Find fa,b,cg maximizing F ¼

Z

qds

ð31Þ

G1 [G2 [G2

Initial design variables are set as {a, b, c}¼{0.10, 0.05, 0.30} m. Area of the fin during iterations, A, is kept equal to the initial area, A0, as an equality constraint as follows: A 1 ¼ 0 A0

ð32Þ

Upper and lower bounds of the design variables are set as follows: 0:02 r a r 0:18 0:02 r b r0:18 0:02 r c r 1:00

ð33Þ

After numerical implementation, the optimal design is obtained in 33 iterations from the initial design and optimal design variables are {a, b, c}¼{0.04420, 0.03853, 0.43454} m. The initial and optimal shapes of the fin are shown in Fig. 13, the

Fig. 14. Temperature distribution of the cooling fin at the initial and optimal designs(unit: 1C). (a) Initial design. (b) Optimal design.

optimal fin being thinner and longer than the initial fin. Values of the objective function F for the initial and optimal designs are 517.6 and 580.4, respectively. Fig. 14(a) and (b) represent temperature distribution for the initial and optimal designs, respectively. It is seen that overall temperature of the wall and fin for the optimal design is lower than the initial design.

5. Conclusions A generalized formulation of the SDSA for two-dimensional steady-state thermal conduction problem as applied to zonedinhomogeneous solids is presented using the BIE formulation and the AVM. The material derivative concept is used to represent the shape variation of the external and interface boundaries. State equation of the primal and adjoint systems is solved by the BEM. Through an analytical example of a composite wall problem, it is proved that the derived sensitivity formula coincides with the analytic solution. For a compound cylinder, a thermal diffuser, and a cooling fin problem, shape sensitivity is numerically

ARTICLE IN PRESS B.Y. Lee / Engineering Analysis with Boundary Elements 34 (2010) 825–833

analyzed using the derived formula and its accuracy is compared with that by the method of finite difference. The sensitivity results appeared fairly accurate overall, but the accuracy of the sensitivity at the elements adjacent to the zone-interface corner is a little lower than elsewhere. As application to numerical optimization, the sensitivity formula is incorporated to a nonlinear programming algorithm and optimum shapes are found for the thermal diffuser problem and the cooling fin problem. Only two-dimensional problems are considered in this work although the formula derivation for axisymmetric or three-dimensional problems can be done in the same way. An implementation study for axisymmetric problems is under way. References [1] Burczynski T. Recent advances in boundary element approach to design sensitivity analysis—a survey. In: Kleiber M, Hisada T, editors. Designsensitivity analysis. Atlanta: Atlanta Technology Publications; 1993. p. 1–25. [2] Lee BY. Direct differentiation formulation for boundary element shape sensitivity analysis of axisymmetric elastic solids. Int J Solids Struct 1997;34: 99–112. [3] Choi JH, Lee BY, Han JS. Boundary integral method for shape optimization of interface problem and its application to implant design in dentistry. Comp Methods Appl Mech Eng 2001;190:6909–26. [4] Li Q, Steven GP, Xie YM, Querin OM. Evolutionary topology optimization for temperature reduction of heat conducting fields. Int J Heat Mass Transfer 2004;47:5071–83. [5] Haug EJ, Choi KK, Komkov V. Design sensitivity analysis of structural systems. New York: Academic Press; 1986. [6] Meric RA. Boundary elements in shape design sensitivity analysis of thermoelastic solids. In: Mota Soares CA, editor. Computer aided optimal design: structural and mechanical systems. Berlin Heidelberg: SpringerVerlag; 1987. p. 7643–52. [7] Meric RA. Shape design sensitivity analysis for nonlinear anisotropic heat conducting solids and shape optimization by the BEM. Int J Numer Methods Eng 1988;26:109–20. [8] Meric RA. Differential and integral sensitivity formulations and shape optimization by BEM. Eng Anal Boundary Elem 1995;15:181–8.

833

[9] Park CW, Yoo YM. Shape design sensitivity analysis of a two-dimensional heat transfer system using the boundary element method. Comp Struct 1988;28:543–50. [10] Lee BY, Kwak BM. Shape optimization of two-dimensional thermoelastic structures using boundary integral equation formulation. Comp Struct 1991;41:709–22. [11] Lee BY, Kwak BM. Axisymmetric thermoelastic shape sensitivity analysis and its application to turbine disc design. Int J Numer Methods Eng 1992;33: 2073–89. [12] Lee BY, Choi JH, Kwak BM. Shape optimization of two-dimensional thermal conducting solid using boundary integral equation formulation. Korean Soc Mech Eng J 1992;6:114–21. [13] Lee BY, Choi JH, Kwak BM. Direct differentiation method for shape design sensitivity analysis of thermal conducting solids. In: Kleiber M, Hisada T, editors. Design-sensitivity analysis. Atlanta: Atlanta Technology Publications; 1993. p. 140–65. [14] Lee BY. Shape sensitivity formulation for axisymmetric thermal conducting solids. Proc Inst Mech Eng, Part C J Mech Eng Sci 1993;207:209–16. [15] Aithal R, Saigal S. Shape sensitivity analysis in thermal problems using BEM. Eng Anal Boundary Elem 1995;15:115–20. [16] Sluzalec A, Kleiber M. Shape sensitivity analysis for nonlinear steady-state heat conduction problems. Int J Heat Mass Transfer 1996;12:2609–13. [17] Kane JH, Wang H. Boundary element shape sensitivity analysis formulations for thermal problems with nonlinear boundary conditions. AIAA J 1991;29: 1978–89. [18] Kane JH, Wang H. Boundary formulations for shape sensitivity of temperature dependent conductivity problems. Int J Numer Methods Eng 1992;33: 2073–2089. [19] Dems K, Mroz Z. Sensitivity analysis and optimal design of external boundaries and interfaces for heat conduction systems. J Thermal Stresses 1998;21:461–88. [20] Gao X-W, He M-C. A new inverse analysis approach for multi-region heat conduction BEM using complex-variable-differentiation method. Eng Anal Boundary Elem 2005;29:788–95. [21] Arora JS. An algorithm for optimum structural design without line search. In: Atrek E, editor. New directions in optimum structural design. New York: A Wiley-Interscience Publications; 1984. p. 429–42. [22] Gao XW, Guo L, Zhang Ch. Three-step multi-domain BEM solver for nonhomogeneous material problems. Eng Anal Boundary Elem 2007;31: 965–73. [23] Mustoe GGW., A combination of the finite element method and boundary integral procedure for continuum problems, Ph.D. thesis, University of Wales, University College, Swansea, 1980.