Special Green’s function boundary element approach for steady-state axisymmetric heat conduction across low and high conducting planar interfaces

Applied Mathematical Modelling 37 (2013) 1948–1965

Contents lists available at SciVerse ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Special Green’s function boundary element approach for steady-state axisymmetric heat conduction across low and high conducting planar interfaces E.L. Chen ⇑, W.T. Ang Division of Engineering Mechanics, School of Mechanical and Aerospace Engineering, Nanyang Technological University, Republic of Singapore

a r t i c l e

i n f o

Article history: Received 16 January 2012 Received in revised form 24 April 2012 Accepted 27 April 2012 Available online 11 May 2012 Keywords: Green’s functions Bimaterial Low and high conducting interfaces Boundary element method

a b s t r a c t The problem of determining the steady-state axisymmetric temperature distribution in a bimaterial with a planar interface is considered here. The interface is either low or high conducting. Special Green’s functions satisfying the thermal conditions on the interface are derived and employed to obtain boundary integral equations whose path of integration does not include the interface. Boundary element procedures that do not require the interface to be discretized into elements are proposed for solving the problem under consideration. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Two dissimilar materials bonded together with a very thin layer of material sandwiched in between them may be modeled as a bimaterial with an interface in the form of a line (for plane problems) or a surface (for three-dimensional problems). For heat conduction problems, the thermal conditions to impose on the line or surface interface depend on the thermal conductivity of the thin layer and may be derived using the asymptotic analysis in [1]. If the temperature and the normal heat flux are continuous on the interface, it (the interface) is considered as ideal (perfect). If the thin (interphase) layer has an extremely low thermal conductivity, the layer may be modeled as an interface over which the temperature is discontinuous. Such an interface is said to be thermally low conducting. For example, the interface between two imperfectly joined solids may be modeled as low conducting if it contains microscopic gaps filled with air. On the other extreme, the thin layer may be modeled as a high conducting interface with discontinuous normal heat flux if the layer is occupied by a material of high thermal conductivity such as carbon nanotubes [2]. Some Green’s functions for steady-state heat conduction across planar interfaces between dissimilar materials may be found in the literature. A Green’s function for two-dimensional heat conduction across an ideal interface between two thermally anisotropic half-spaces was used in [3] to develop a meshless method for analyzing the temperature distribution in a bimaterial. Ang et al. [4] derived a Green’s function for two-dimensional heat conduction across a low conducting interface between two thermally isotropic half-spaces. The Green’s function was applied to derive boundary element procedures for solving two-dimensional heat conduction problems involving bimaterials of finite extent (see also [5]). Three-dimensional Green’s functions for heat conduction across low and high conducting interfaces were given in [6]. ⇑ Corresponding author. E-mail address: [email protected] (E.L. Chen). 0307-904X/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.apm.2012.04.051

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

1949

In the present paper, we consider the problem of determining the steady-state axisymmetric temperature distribution in an axisymmetric bimaterial with a planar interface which is either low or high conducting. For thermal analysis of the bimaterial, special Green’s functions satisfying the relevant thermal conditions on the interface are derived and employed to obtain boundary integral equations whose path of integration does not include the interface. Boundary element procedures that do not require the interface to be discretized into elements are proposed. To check the validity and accuracy of the boundary element procedures, they are applied to solve some specific problems. 2. The problem Consider two dissimilar materials bonded together with a thin layer of material sandwiched in between them. The regions occupied by the layer and the two dissimilar materials are denoted by R0 ; R1 and R2 , respectively. With reference to a Cartesian coordinate system denoted by Oxyz; R0 occupies part of the space h=2 < z < h=2, where h is a given positive number. The regions R1 and R2 are subsets of the half-spaces z < h=2 and z > h=2, respectively. The regions R0 ; R1 and R2 are axisymmetric, obtained by rotating respectively two-dimensional regions X0 ; X1 and X2 on the Orz (axisymmetric coordinate) plane (as sketched in Fig. 1) by an angle of 360 about the z-axis. Note that r is the distance of a point from the z-axis. The thermal conductivities of the materials in X0 ; X1 and X2 are positive constants j0 ; j1 and j2 , respectively. We are interested in modeling the sandwiched layer X0 as a line interface on the z-axis of the Orz plane for the limiting case in which the thickness h tends to zero. A geometrical sketch of the body with the line interface a < r < b; z ¼ 0 denoted by C0 is shown in Fig. 2. (Note that the constant a is 0 if X0 and C0 are as sketched in Figs. 1 and 2, respectively. In general, depending on the geometries of X0 ; X1 and X2 ; a is not necessarily zero.) The asymptotic approach in [1] may be employed to derive appropriate thermal conditions to impose on the line interface C0 . It is assumed here that the temperature in the materials is steady and varies spatially with the axisymmetric coordinates r and z only. Denoted by TðxÞ (where x ¼ ðr; zÞ), the temperature satisfies the governing partial differential equation

@ 2 T 1 @T @ 2 T þ þ ¼ 0 for x 2 X0 [ X1 [ X2 : @r 2 r @r @z2

ð1Þ

The line interface C0 between X1 and X2 in Fig. 2 is said to be ideal or perfectly conducting if the temperature T and the normal heat flux are continuous on C0 , that is, if the thermal conditions on the interface are given by

Tðr; 0þ Þ ¼ Tðr; 0 Þ   j2 @T  þ ¼ j1 @T  @z z¼0

) for a < r < b:

ð2Þ

@z z¼0

If the thermal conductivity

j0 in the layer X0 is such that

Fig. 1. The regions X0 ; X1 and X2 are rotated by an angle of 360 about the z-axis to form the three-dimensional regions R0 ; R1 and R2 , respectively.

1950

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

Fig. 2. The layer X0 is replaced by the line C0 as h tends to zero. Apart from C0 , the curves C1 and C2 make up the remaining boundaries of X1 and X2 , respectively.

j0 h

! k ða finite positive constantÞ as h ! 0þ ;

ð3Þ

then the asymptotic analysis in [1] may be used to derive the following thermal conditions on C0 :



)



  j2 @T ¼ j1 @T @z z¼0þ @z z¼0

 k½Tðr; 0þ Þ  Tðr; 0 Þ ¼ j2 @T 

for a < r < b:

ð4Þ

@z z¼0þ

Note that (3) implies that j0 approaches zero as h tends to zero. Thus, (4) gives the thermal conditions on a layer with thickness and thermal conductivity that tend to zero in such a way that there is a temperature jump across opposite sides of the layer of vanishing thickness. A non-ideal conducting interface with such thermal conditions is said to be low conducting. The line interface C0 between two materials as sketched in Fig. 2 may be modeled as low conducting if the interface is imperfect containing microscopic gaps filled with air. For another non-ideal conducting interface, if the thermal conductivity of the layer X0 is given by

j0 h ! a ða finite positive constantÞ as h ! 0þ ;

ð5Þ

then the thermal conditions on C0 are given by

Tðr; 0þ Þ ¼ Tðr; 0 Þ   j2 @T  þ  j1 @T 

@z z¼0

@z z¼0

 2  ¼ a@@z2T 

z¼0

9 = ;

for a < r < b:

ð6Þ

In view of the first condition in (6) and the governing partial differential equation in (1), note that

 @ 2 T   @z2 

¼ z¼0þ

     1 @ @T  1 @ @T  @ 2 T  r r ¼  ¼  r @r @r z¼0þ r @r @r z¼0 @z2 

:

ð7Þ

z¼0

As implied by (5), the thermal conductivity j0 in the layer X0 tends to infinity as h vanishes. Thus, a non-ideal conducting interface with the thermal conditions (6) in which the normal heat flux is discontinuous across the vanishing interphase layer is said to be high conducting. For a practical example, if the two materials in Fig. 2 are joined together by an extremely thin layer of carbon nanotubes, the line interface may be modeled as high conducting. The problem of interest here is to solve (1) for the axisymmetric steady-state temperature in the bimaterial sketched in Fig. 2, that is, in X1 [ X2 , subject to either (4) or (6) (as the thermal conditions on C0 ) and suitably prescribed temperature or flux at each point on the exterior boundary C1 [ C2 of the bimaterial. Specifically, the boundary conditions on C1 [ C2 are given by

TðxÞ ¼ f0 ðxÞ for x 2 N1 ; Pðx; nðxÞÞ ¼ f1 ðxÞ þ f2 ðxÞTðxÞ for x 2 N2 ;

ð8Þ

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

1951

where f0 ðxÞ; f1 ðxÞ and f2 ðxÞ are suitably prescribed functions, N1 and N2 are non-intersecting curves (for the different boundary conditions) such that N1 [ N2 ¼ C1 [ C2 . 3. Green’s function boundary element method 3.1. Boundary integral equations The boundary integral equations for axisymmetric heat conduction governed by (1) in the bimaterial sketched in Fig. 2 are given by (see [7]) Z

c1 ðx0 ÞTðx0 Þ ¼

ðTðxÞG1 ðx; x0 ; nðxÞÞ  G0 ðx; x0 ÞPðx; nðxÞÞÞrdsðxÞ

C1

Z

þ

b a

and

c2 ðx0 ÞTðx0 Þ ¼

    @ rdr Tðr; 0 ÞG1 ðr; 0 ; x0 ; 0; 1Þ  G0 ðr; 0 ; x0 Þ ½TðxÞ @z z¼0

for x0 ¼ ðr0 ; z0 Þ 2 X1 [ C0 [ C1

ð9Þ

Z

ðTðxÞG1 ðx; x0 ; nðxÞÞ  G0 ðx; x0 ÞPðx; nðxÞÞÞrdsðxÞ   Z b  @  rdr Tðr; 0þ ÞG1 ðr; 0þ ; x0 ; 0; 1Þ  G0 ðr; 0þ ; x0 Þ ½TðxÞ @z a z¼0þ C2

for x0 2 X2 [ C0 [ C2 ;

ð10Þ

where ci ðx0 Þ ¼ 1 if x0 lies in the interior of Xi ; c1 ðx0 Þ and c2 ðx0 Þ are defined by

c1 ðx0 Þ ¼ and

c2 ðx0 Þ ¼

Z

G1 ðx; x0 ; nðxÞÞrdsðxÞ þ

b

G1 ðr; 0 ; x0 ; 0; 1Þrdr

for x0 ¼ ðr 0 ; z0 Þ 2 X1 [ C0 [ C1

ð11Þ

G1 ðr; 0þ ; x0 ; 0; 1Þrdr

for x0 2 X2 [ C0 [ C2 ;

ð12Þ

a

C1

Z

Z

G1 ðx; x0 ; nðxÞÞrdsðxÞ 

Z

b

a

C2

dsðxÞ denotes the length of an infinitesimal part of the curve C0 [ Ci ; nðxÞ ¼ ½nr ðxÞ; nz ðxÞ ¼ nr ðxÞer þ nz ðxÞez (er and ez are the unit base vectors along the r and z axes, respectively) is the unit normal vector to C1 [ C2 (at the point x) pointing out of X1 [ X2 ; Pðx; nðxÞ) is the directional rate of change of the axisymmetric temperature along the vector n(x) as defined by

Pðx; nðxÞÞ ¼ nr ðxÞ

@ @ ½TðxÞ þ nz ðxÞ ½TðxÞ @r @z

ð13Þ

and G0 ðx; x0 Þ and G1 ðx; x0 ; nðxÞÞ are given by ð1Þ

ð2Þ

G0 ðx; x0 Þ ¼ ½HðzÞHðz0 Þ þ HðzÞHðz0 ÞG0 ðx; x0 Þ þ G0 ðx; x0 Þ; ð1Þ

ð2Þ

G1 ðx; x0 ; nðxÞÞ ¼ ½HðzÞHðz0 Þ þ HðzÞHðz0 ÞG1 ðx; x0 ; nðxÞÞ þ G1 ðx; x0 ; nðxÞÞ; Kðmðx; x0 ÞÞ ð1Þ G0 ðx; x0 Þ ¼  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; p aðx; x0 Þ þ bðr; r0 Þ 1 ð1Þ G1 ðx; x0 ; nðxÞÞ ¼  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p aðx; x0 Þ þ bðr; r0 Þ ( " # ) nr ðxÞ r 20  r 2 þ ðz0  zÞ2 z0  z Eðmðx; x0 ÞÞ  Kðmðx; x0 ÞÞ þ nz ðxÞ Eðmðx; x0 ÞÞ ;  2r aðx; x0 Þ  bðr; r0 Þ aðx; x0 Þ  bðr; r0 Þ ð14Þ with H denoting the unit-step Heaviside function, the functions mðx; x0 Þ; aðx; x0 Þ; KðmÞ and EðmÞ given by

2bðr; r 0 Þ ; aðx; x0 Þ þ bðr; r 0 Þ 2 2 2 aðx; x0 Þ ¼ r 0 þ r þ ðz0  zÞ ; bðr; r 0 Þ ¼ 2rr 0 ; Z p=2 dh pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; KðmÞ ¼ 2 0 1  m sin h Z p=2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 EðmÞ ¼ 1  m sin hdh; mðx; x0 Þ ¼

0 ð2Þ

the function G0 ðx; x0 Þ being any solution of

ð15Þ

1952

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

i 1 @ h i @2 h i @ 2 h ð2Þ ð2Þ ð2Þ G0 ðx; x0 Þ þ G0 ðx; x0 Þ þ 2 G0 ðx; x0 Þ ¼ 0 for x 2 X1 [ X2 2 r @r @r @z

ð16Þ

ð2Þ

and the function G1 ðx; x0 ; nðxÞÞ defined by ð2Þ

G1 ðx; x0 ; nðxÞÞ ¼ nr ðxÞ We may multiply

i i @ h ð2Þ @ h ð2Þ G0 ðx; x0 Þ þ nz ðxÞ G0 ðx; x0 Þ : @r @z

ð17Þ

j1 and j2 to (9) and (10), respectively and add up the two equations to obtain

c1 ðx0 Þj1 Tðx0 Þ þ c2 ðx0 Þj2 Tðx0 Þ ¼

2 Z X i¼1

þ

Z

ji ðTðxÞG1 ðx; x0 ; nðxÞÞ  G0 ðx; x0 ÞPðx; nðxÞÞÞrdsðxÞ

Ci b

a



Z

a

b



 

@ @z

j1 Tðr; 0 ÞG1 ðr; 0 ; x0 ; 0; 1Þ  G0 ðr; 0 ; x0 Þ ½TðxÞ

 rdr

z¼0

    @ rdr j2 Tðr; 0þ ÞG1 ðr; 0þ ; x0 ; 0; 1Þ  G0 ðr; 0þ ; x0 Þ ½TðxÞ @z z¼0þ

for x0 2 X1 [ X2 [ C1 [ C2 :

ð18Þ

ð2Þ G0 ðx; x0 Þ

ð2Þ G0 ðx; x0 Þ

In general, we may take in (14) to be ¼ 0. Nevertheless, we may find it advantageous to solve (16) for ð2Þ G0 ðx; x0 Þ that satisfies certain thermal conditions on the interface C0 of the bimaterial. As we shall show below, a specially ð2Þ chosen G0 ðx; x0 Þ may be used in (14) for the integral Eq. (18) such that the integrals over the interface C0 vanishes. 3.2. Green’s function for low conducting interface ð2Þ

For the case in which the interface C0 of the bimaterial is low conducting, G0 ðx; x0 Þ is chosen in such a way that G0 ðx; x0 Þ in (14) satisfies the interfacial conditions



)



j2 @G@z0 z¼0þ ¼ j1 @G@z0 z¼0

 0 k½G0 ðr; 0þ ; x0 Þ  G0 ðr; 0 ; x0 Þ ¼ j2 @G @z z¼0þ

for 0 < r < 1:

ð19Þ

The Green’s function G0 ðx; x0 Þ satisfying (19) can be obtained by performing an axial integration on the corresponding Green’s function for three-dimensional heat conduction across a low conducting planar interface at z = 0. The analysis in [6] may be easily adapted to derive the corresponding three-dimensional Green’s function. The derivation is given in Appendix A. If we let x ¼ r cos h; y ¼ r sin h; n ¼ r 0 ; g ¼ 0 and f ¼ z0 in the three-dimensional Green’s function for the low conducting interface and integrate it with respect to h from h = 0 to h ¼ 2p, we find that the required axisymmetric Green’s function G0 ðx; x0 Þ for low conducting C0 is given by (14) with

  Z 2k ð2Þ ð1Þ G0 ðx; x0 Þ ¼ Hðz0 Þ HðzÞ G0 ðx; r 0 ; z0 Þ  þ HðzÞ

2k

j1

Z

j2

1

0

þ HðzÞ

Z

j1

0

     k j2 ð1Þ u du G0 ðx; r 0 ; u  z0 Þ exp  1þ

    k j2 ð1Þ u du G0 ðx; r 0 ; z0  uÞ exp  1þ

  Z 2k ð1Þ þ Hðz0 Þ HðzÞ G0 ðx; r 0 ; z0 Þ  2k

1

j2

0

j2

1

     k j2 ð1Þ u du G0 ðx; r 0 ; z0  uÞ exp  1þ

j2

0

j1

j1

    k j2 ð1Þ u du : G0 ðx; r 0 ; z0 þ uÞ exp  1þ

1

j2

j2

j1

ð20Þ

j1

ð2Þ

Note that G0 ðx; x0 Þ is obtained by integrating axially U ðx; y; z; n; g; fÞ given in Appendix A by (A4) and (A11). ð2Þ ð2Þ The function G1 ðx; x0 ; nðxÞÞ which corresponds to G0 ðx; x0 Þ in (20) is given by

  Z 2k ð2Þ ð1Þ G1 ðx; x0 ; nðxÞÞ ¼ Hðz0 Þ HðzÞ G1 ðx; r 0 ; z0 ; nðxÞÞ  þHðzÞ

2k

j1

Z

j2

1

0

0



þHðzÞ

2k

j1

0

   k j2 ð1Þ G1 ðx; r 0 ; u  z0 ; nðxÞÞ exp  ð1 þ Þu du

 k j2 ð1Þ G1 ðx; r 0 ; z0  u; nðxÞÞ exp  ð1 þ Þu du

  Z 2k ð1Þ þ Hðz0 Þ HðzÞ G1 ðx; r 0 ; z0 ; nðxÞÞ  Z

1

j2

1

0

j2

1

j2

j1

j1

     k j2 ð1Þ u du G1 ðx; r0 ; z0  u; nðxÞÞ exp  1þ

  k j2 ð1Þ G1 ðx; r 0 ; z0 þ u; nðxÞÞ exp  ð1 þ Þu du :

j2

j1

Using (4) and (19) for low conducting C0 , we find that (18) can be reduced to

j2

j1

ð21Þ

1953

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

fc1 ðx0 Þj1 þ c2 ðx0 Þj2 gTðx0 Þ ¼

2 Z X

ji ðTðxÞG1 ðx; x0 ; nðxÞÞ  G0 ðx; x0 ÞPðx; nðxÞÞÞrdsðxÞ for x0 2 X1 [ X2 [ C1 [ C2 ;

Ci

i¼1

ð22Þ if G0 ðx; x0 Þ and G1 ðx; x0 ; nðxÞÞ are given by (14) with (20) and (21). If (22) is used to derive a boundary element procedure for the numerical solution of the problem stated in Section 2, it is not necessary to discretize the low conducting interface C0 into boundary elements. Thus, the system of linear algebraic equations in the boundary element formulation is smaller with fewer unknowns. 3.3. Green’s function for high conducting interface ð2Þ

For the case in which the interface C0 of the bimaterial is high conducting, G0 ðx; x0 Þ is chosen in such a way that G0 ðx; x0 Þ in (14) satisfies the interfacial conditions

9 =

G0 ðr; 0þ ; x0 Þ ¼ G0 ðr; 0 ; x0 Þ    2 j2 @G@z0 z¼0þ  j1 @G@z0 z¼0 ¼ a@@zG20 

z¼0

for 0 < r < 1:

;

ð23Þ

ð2Þ

The function G0 ðx; x0 Þ such that (23) holds is obtained by integrating axially U ðx; y; z; n; g; fÞ given by (A4) and (A19) in Appendix A (for three-dimensional heat conduction across a high conducting planar interface at z = 0), that is,

  Z 2j1 ð2Þ ð1Þ G0 ðx; x0 Þ ¼ Hðz0 Þ HðzÞ G0 ðx; r 0 ; z0 Þ þ

a

þ HðzÞ

2j1

Z

a

0

0

a

a

a

Z

2j2

   1 ð1Þ G0 ðx; r 0 ; u  z0 Þ exp  ðj1 þ j2 Þu du

  1 1 ð1Þ G0 ðx; r0 ; z0  uÞ exp  ðj1 þ j2 Þu du

  Z 2j2 ð1Þ þ Hðz0 Þ HðzÞ G0 ðx; r 0 ; z0 Þ þ

þ HðzÞ

1

a

1

1

0

   1 ð1Þ G0 ðx; r 0 ; z0  uÞ exp  ðj1 þ j2 Þu du

a

  1 ð1Þ G0 ðx; r 0 ; z0 þ uÞ exp  ðj1 þ j2 Þu du :

ð24Þ

a

0

ð2Þ

ð2Þ

The function G1 ðx; x0 ; nðxÞÞ which corresponds to G0 ðx; x0 Þ in (24) is given by

  Z 2j1 ð2Þ ð1Þ G1 ðx; x0 ; nðxÞÞ ¼ Hðz0 Þ HðzÞ G1 ðx; r0 ; z0 ; nðxÞÞ þ

a

þHðzÞ

2j1

Z

a

0

a

0

a

a

2j2

   1 ð1Þ G1 ðx; r0 ; u  z0 ; nðxÞÞ exp  ðj1 þ j2 Þu du

  1 1 ð1Þ G1 ðx; r 0 ; z0  u; nðxÞÞ exp  ðj1 þ j2 Þu du

  Z 2j2 ð1Þ þ Hðz0 Þ HðzÞ G1 ðx; r0 ; z0 ; nðxÞÞ þ

þHðzÞ

1

Z

a

0

0

1

 1 ð1Þ G1 ðx; r 0 ; z0  u; nðxÞÞ expð ðj1 þ j2 ÞuÞdu

a

  1 1 ð1Þ G1 ðx; r 0 ; z0 þ u; nðxÞÞ exp  ðj1 þ j2 Þu du :

a

ð25Þ

For high conducting C0 , using (6), (7) and (23) and noting that

 @ 2 G0   @z2 

¼ z¼0

  1 @ @ r ½G0 ðr; 0; x0 Þ ; r @r @r

ð26Þ

we find that (18) can be rewritten as

c1 ðx0 Þj1 Tðx0 Þ þ c2 ðx0 Þj2 Tðx0 Þ ¼

2 Z X

þ

ji ðTðxÞG1 ðx; x0 ; nðxÞÞ  G0 ðx; x0 ÞPðx; nðxÞÞÞrdsðxÞ

Ci

i¼1

Z a

for x0 2 X1 [ X2 [ C1 [ C2 : Using integration by parts, we find that

b



a G0 ðr; 0; x0 Þ

    1 @ @T  1 @ @ r r ½G þ Tðr; 0Þ ðr; 0; x Þ rdr 0 0 r @r @r z¼0þ r @r @r ð27Þ

1954

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

    1 @ @T  1 @ @ r r ½G þ Tðr; 0Þ ðr; 0; x Þ rdr 0 0 r @r @r z¼0þ r @r @r a     @T  @T  @ þ aG0 ða; 0; x0 Þa  þ bTðb; 0Þa ½G0 ðx; x0 Þ ¼ bG0 ðb; 0; x0 Þa  @r ðr;zÞ¼ðb;0Þ @r ðr;zÞ¼ða;0Þ @r ðr;zÞ¼ðb;0Þ   @  aTða; 0Þa ½G0 ðx; x0 Þ : @r ðr;zÞ¼ða;0Þ

Z

b



a G0 ðr; 0; x0 Þ

ð28Þ

It follows that (27) reduces to

    @ @ þ aTða; 0Þa ½G0 ðx; x0 Þ fc1 ðx0 Þj1 þ c2 ðx0 Þj2 gTðx0 Þ  bTðb; 0Þa ½G0 ðx; x0 Þ @r @r ðr;zÞ¼ðb;0Þ ðr;zÞ¼ða;0Þ   @T  @T  þ bG0 ðb; 0; x0 Þa   aG0 ða; 0; x0 Þa  @r ðr;zÞ¼ðb;0Þ @r ðr;zÞ¼ða;0Þ Z 2 X ¼ ji ðTðxÞG1 ðx; x0 ; nðxÞÞ  G0 ðx; x0 ÞPðx; nðxÞÞÞrdsðxÞ for x0 2 X1 [ X2 [ C1 [ C2 : i¼1

ð29Þ

Ci

Thus, the integral over the high conducting interface C0 vanishes if the Green’s function G0 ðx; x0 Þ is given by (14) with ð2Þ G0 ðx; x0 Þ in (25). 3.4. Boundary element procedures In this section, we describe boundary element procedures for determining TðxÞ and Pðx; nðxÞÞ (whichever is not known) on C1 [ C2 . Once TðxÞ and Pðx; nðxÞÞ are completely known on C1 [ C2 , we can obtain the temperature at any point x0 in the interior of the domains by using c1 ðx0 Þ ¼ 1 and c2 ðx0 Þ ¼ 0 for x0 in the interior of X1 or c1 ðx0 Þ ¼ 0 and c2 ðx0 Þ ¼ 1 for x0 in the interior of X2 in (22) (for low conducting interface) or in (29) (for high conducting interface). We discretize the boundary C1 [ C2 into N straight line elements denoted by Bð1Þ ; Bð2Þ ;    ; BðN1Þ and BðNÞ . As (22) or (29) does not contain any integral over the low conducting interface (because of the use of the special Green’s function), we do not need to discretize the interface C0 . For a simple approximation, T and P are taken to be constants over an element of C1 [ C2 , specifically

TðxÞ ’ T ðmÞ

) for x 2 BðmÞ

Pðx; nðxÞÞ ’ PðmÞ

ðm ¼ 1; 2; . . . ; NÞ;

ð30Þ

where T ðmÞ and P ðmÞ are constants. Each boundary element is associated with only one unknown constant. Specifically, if T is specified over the element BðmÞ according to the first line of (8) then P ðmÞ is the unknown over BðmÞ . On the other hand, if P is given by the second line of (8) over BðmÞ , we can express P ðmÞ in terms of T ðmÞ and regard T ðmÞ as the unknown constant over BðmÞ . 3.4.1. Low conducting interface For low conducting interface C0 , letting x0 in (22) be given in turn by the midpoints of BðiÞ ði ¼ 1; 2; . . . ; NÞ, together with (8), we obtain

n

ðiÞ

ðiÞ

oh

ðiÞ

i

c1 ð xb Þj1 þ c2 ð bx Þj2 dðiÞ T ðiÞ þ ð1  dðiÞ Þf0 ð xb Þ ¼

N X

nh

jðmÞ dðmÞ T ðmÞ þ ð1  dðmÞ Þf0 ð xb ðmÞ Þ

i

m¼1



Z BðmÞ



h

 ðiÞ ðmÞ b ðmÞ Þ þ f2 ð x b ðmÞ ÞT ðmÞ x ; nðmÞ ÞrdsðxÞ  d G1 ðx; b f1 ð x ðmÞ

þ 1d ðiÞ

ðmÞ

ÞPðmÞ

iZ

ðiÞ

BðmÞ

b ÞrdsðxÞ G0 ðx; x

for i ¼ 1; 2; . . . ; N; ðmÞ

ð31Þ

b is the midpoint of BðiÞ ; d ¼ 0 if T is specified on the mth element BðmÞ as given by the first line of (8), d ¼ 1 if the where x boundary condition given by the second line of (8) is applicable on BðmÞ ; nðmÞ is the unit normal vector to BðmÞ pointing away from the solution domain X1 [ X2 ; jðmÞ ¼ j1 if BðmÞ is an element on the boundary of X1 and jðmÞ ¼ j2 if BðmÞ is an element on the boundary of X2 . b ðiÞ is the midpoint of BðmÞ (that is, if m = i). The Cauchy principal In (31), the integrals over BðmÞ are Cauchy principal if x integrals can be accurately evaluated by using a highly accurate Gaussian quadrature. Now (31) gives a system of N linear algebraic equations containing N unknowns given by either T ðkÞ or P ðkÞ ðk ¼ 1; 2; . . . ; NÞ. Once the unknowns on the boundary are determined, the temperature at the interior point of the domain X1 [ X2 can be obtained as explained above.

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

1955

3.4.2. High conducting interface For high conducting interface C0 , if we proceed as before by collocating (29) at the midpoint of each boundary element, we obtain

n

ðiÞ

oh

ðiÞ

ðiÞ

i

@ @r

 

c1 ð xb Þj1 þ c2 ð xb Þj2 dðiÞ T ðiÞ þ ð1  dðiÞ Þf0 ð xb Þ  bTðb; 0Þa ½G0 ðx; x0 Þ

ðr;zÞ¼ðb;0Þ

    @ @T  @T  þ aTða; 0Þa ½G0 ðx; x0 Þ þ bG0 ðb; 0; x0 Þa   aG0 ða; 0; x0 Þa  @r @r ðr;zÞ¼ðb;0Þ @r ðr;zÞ¼ða;0Þ ðr;zÞ¼ða;0Þ  Z N X b ðiÞ ; nðmÞ ÞrdsðxÞ  ½dðmÞ ðf1 ð x b ðmÞ Þ ¼ jðmÞ ½dðmÞ T ðmÞ þ ð1  dðmÞ Þf0 ð bx ðmÞ Þ G1 ðx; x m¼1

b ðmÞ ÞT ðmÞ Þ þ ð1  dðmÞ ÞP ðmÞ  þf2 ð x

Z BðmÞ

BðmÞ

b ðiÞ ÞrdsðxÞ G0 ðx; x

for i ¼ 1; 2; . . . ; N;

ð32Þ

where jðmÞ is as defined below  (31).    The terms Tða; 0Þ; Tðb; 0Þ; @T and @T in (32) are unknown constants. They can, however, be approximated in @r ðr;zÞ¼ða;0Þ @r ðr;zÞ¼ðb;0Þ terms of T and P on boundary elements near (a, 0) and/or (b, 0). How the required approximations may be made depends on the geometries of the solution domains–see, for example, Problems 3 and 4 in Section 4 below. Thus, (32) can be solved as a system of N linear algebraic equations for N unknowns given by either T ðkÞ or P ðkÞ ðk ¼ 1; 2; . . . ; NÞ. 4. Specific problems

Problem 1. To test the boundary element procedure for C0 that is low conducting, consider the regions X1 and X2 as sketched in Fig. 3. Note that X1 and X2 are defined by the curves r2 þ z2 ¼ 4 and r 2 þ z2 ¼ 1 and the lines r ¼ 0; r ¼ 1; r ¼ 2 and z ¼ 1 on the rz plane. For a particular problem take j1 ¼ 1; j2 ¼ 2 and k ¼ 1. The exterior boundary of X1 [ X2 is approximated using N straight line elements. The boundary conditions on the exterior boundary of X1 [ X2 are given by

Pð2; z; 1; 0Þ ¼ 4z Pð1; z; 1; 0Þ ¼ 2z



for

 1 < z < 0;

2 for 1 < r < 2; 3 1 1 Tðr; zÞ ¼ r 2 z  z3 þ r 2  2z2 for r 2 þ z2 ¼ 1; 0 < r < 1; 2  3  1 1 1 3 P r; z; r; z ¼  z3  2z2 þ r2 z þ r 2 for r 2 þ z2 ¼ 4; 0 < r < 2: 2 2 2 4 Tðr; 1Þ ¼ r 2 þ

Fig. 3. A geometrical sketch of Problem 1 on the rz plane.

1956

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

It can be easily verified that the exact solution for the problem here is given by

( Tðr; zÞ ¼

r 2 z  23 z3  13 z3 þ r 2  2z2

1 2 r z 2

for ðr; zÞ 2 X1 ; for ðr; zÞ 2 X2 :

Numerical values are obtained for the temperature at the interior of the regions by solving Eq. (31) using N = 60 and N = 120. Through the use of (22) with c1 ðx0 Þ ¼ 1; c2 ðx0 Þ ¼ 0 for x0 in the interior of X1 and c1 ðx0 Þ ¼ 0; c2 ðx0 Þ ¼ 1 for x0 in the interior of X2 , numerical values of the temperature at r = 1.5 and r = 1.75 for 1 < z < 0 and at r2 þ z2 ¼ ð1:75Þ2 and r2 þ z2 ¼ ð1:5Þ2 for 0 < z < 1:5 are obtained and compared graphically with the exact temperature in Fig. 4. On the whole, the numerical and exact temperature agree well with each other. Note that the gap in the graph is due to the temperature jump across the interface C0 at z = 0. Problem 2. Consider now the case in which X1 and X2 are given by

X1 ¼ fðr; zÞ : 0 6 r < 1; 1 < z < 0g; 

3 2

X2 ¼ ðr; zÞ : 0 6 r < ; 0 < z <

3 ; 2

as sketched in Fig. 5. As in Problem 1, the interface C0 between X1 and X2 is taken to be low conducting. Note that for this particular case the exterior boundary of the bimaterial lies on part of the z = 0 plane (that is, 1 < r < 3=2; z ¼ 0). For a particular problem, we take j1 ¼ 1; j2 ¼ 1=2 and k = 1 and the boundary conditions as

Tð1; zÞ ¼ 4  2z2 þ 2z for

 1 < z < 0;

Pðr; 1; 0; 1Þ ¼ 6 for 0 < r < 1; 3 13 3 for 0 < r < ; Tðr; Þ ¼ r 2 þ 2 2 2 3 3 Pð ; z; 1; 0Þ ¼ 3 for 0 < z < ; 2 2 3 2 Tðr; 0Þ ¼ r þ 5 for 1 < r < : 2 The exact solution for the particular problem here is given by

( Tðr; zÞ ¼

r 2  2z2 þ 2z þ 3 for ðr; zÞ 2 X1 ; r 2  2z2 þ 4z þ 5 for ðr; zÞ 2 X2 :

Fig. 4. Plots of the numerical and exact temperature Tðr; zÞ for 1 < z < 1:5 (Problem 1).

1957

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

Fig. 5. A geometrical sketch of Problem 2 on the rz plane.

The exterior boundary of the bimaterial is discretized into N straight line elements. The numerical values of T at various selected points in X1 [ X2 are computed using (22) with c1 ðx0 Þ ¼ 1 and c2 ðx0 Þ ¼ 0 for x0 in the interior of X1 and c1 ðx0 Þ ¼ 0 and c2 ðx0 Þ ¼ 1 for x0 in the interior of X2 . They are compared with the exact values in Table 1 for N = 25, 50, 100 and 200. The numerical values are reasonably accurate and they converge to the exact solution when the calculation is refined by reducing the sizes of the boundary elements used (that is, when N is increased from 25 to 200). All percentage errors of the numerical values for N = 200 are less than 0.05%. Problem 3. To check the boundary element procedure for a bimaterial with a high conducting interface, take

X1 ¼ fðr; zÞ : 0 6 r < 1; 1 < z < 0g; X2 ¼ fðr; zÞ : 0 6 r < 1; 0 < z < 1g as illustrated in Fig. 6. We take j1 ¼ 6; j2 ¼ 2, and a ¼ 7=4 and the boundary conditions as

Pð1; z; 1; 0Þ ¼ 4 þ z for  1 < z < 0; 1 Tð1; zÞ ¼ 2 þ z  4z2  z3 for 0 < z < 1; 2 ) Tðr; 1Þ ¼ 32 r2  17 3 for 0 < r < 1: Pðr; 1; 0; 1Þ ¼ 32 r 2  12 The exact solution of the particular problem here is given by

( Tðr; zÞ ¼

2r 2  4z2 þ 12 r 2 z  13 z3 þ 2z for ðr; zÞ 2 X1 ; 2r 2  4z2 þ 32 r 2 z  z3  z

for ðr; zÞ 2 X2 :

Table 1 Numerical and exact values of T at selected interior points for Problem 2. Point

N = 25

N = 50

N = 100

N = 200

Exact

(0.30, 0.80) (0.75, 0.35) (0.40, 0.40) (1.40, 0.20) (0.50, 1.25) (0.20, 0.80)

0.21635 2.61420 2.03921 7.60894 7.11813 6.94401

0.21094 2.61655 2.03947 7.66118 7.12282 6.95498

0.21005 2.61724 2.03977 7.67354 7.12431 6.95848

0.20996 2.61743 2.03992 7.67800 7.12478 6.95954

0.21000 2.61750 2.04000 7.68000 7.12500 6.96000

1958

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

Fig. 6. A geometrical sketch of Problem 3 on the rz plane.

The exterior boundary of the bimaterial is discretized into N straight line  elements. To solve (32) as a system of N linear  algebraic equations in N unknowns, we have to approximate T (1, 0) and @T in terms of T and P on the boundary ele@r ðr;zÞ¼ð1;0Þ ments. (Note that for this particular problem, a = 0 and b = 1.) If the exterior boundary is discretized in such a way that the first and the last elements (Bð1Þ and BðNÞ , respectively) are of equal length, lie on r = 1 and have (1, 0) as one of their endpoints, then we can make the approximations

1 Tð1; 0Þ ’ ðT ð1Þ þ T ðNÞ Þ; 2  @T  1 ’ ðPð1Þ þ PðNÞ Þ: @r ðr;zÞ¼ð1;0Þ 2 Numerical values of the temperature at selected interior points are compared with the exact values in Table 2. On the whole, the numerical values at the selected interior points are in good agreement with the exact solution and improve in accuracy when N is increased from 40 to 320 (again, calculation is refined by reducing the sizes of the boundary elements   used). The numerical results here also justify the above approximations for the terms Tð1; 0Þ and @T in (32). @r ðr;zÞ¼ð1;0Þ Problem 4. Problem 3 deals with relatively simple rectangular domains on the rz-plane. For a more general test problem involving a high conducting interface, we take here X1 with a slanted boundary r ¼ z þ 3=2 and X2 with part of z = 0 as its exterior boundary. A sketch of X1 [ X2 is given in Fig. 7. We take j1 ¼ 1; j2 ¼ 1=2 and a ¼ 1=8 and the boundary conditions on the exterior boundary of X1 [ X2 as

1 1 3 Tðr; zÞ ¼ r2  2z2 þ r2 z  z3 for r ¼ z þ ; 1 < z < 0; 2 3 2 1 2 1 Pðr; 1; 0; 1Þ ¼ 3  r for 0 < r < ; 2 2 Pðr; 2; 0; 1Þ ¼ 17 þ r 2 for 0 < r < 2; Pð2; z; 1; 0Þ ¼ 4ð1 þ zÞ for 0 < z < 2; 3 < r < 2: Tðr; 0Þ ¼ r 2 for 2 Table 2 Numerical and exact values of T at selected interior points for Problem 3. Point

N = 40

N = 80

N = 160

N = 320

Exact

(0.500, 0.500) (0.200, 0.100) (0.700, 0.950) (0.900, 0.900) (0.100, 0.500) (0.750, 0.001)

1.52954 0.17288 4.47882 2.14050 1.59994 1.10860

1.52386 0.16598 4.47771 2.15158 1.59895 1.11827

1.52198 0.16350 4.47724 2.15451 1.59828 1.12191

1.52130 0.16249 4.47707 2.15525 1.59790 1.12347

1.52083 0.16167 4.47696 2.15550 1.59750 1.12484

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

1959

Fig. 7. A geometrical sketch of Problem 4 on the rz plane.

The exact solution of the problem here is given by

( Tðr; zÞ ¼

r 2  2z2 þ 12 r 2 z  13 z3 2

2

2

r  2z þ r z 

2 3 z 3

for ðr; zÞ 2 X1 ;

 z for ðr; zÞ 2 X2 :

To solve (32), we have to approximate Tð3=2; 0Þ and



@T  @r ðr;zÞ¼ð3=2;0Þ

in terms of T and P on the boundary elements which

approximate the exterior boundary of X1 [ X2 . If there is a small horizontal boundary element at the intersection between   the slanted boundary r ¼ z þ 3=2 and the vertical boundary z ¼ 0 then @T can be approximated as P on that horizon@r ðr;zÞ¼ð3=2;0Þ tal element. For this purpose, we approximate the line segment r ¼ z þ 3=2; d < z < 0, where d is a very small number, by using a small horizontal line element of length d (see Fig. 7). If the small horizontal line element is taken to be the first element Bð1Þ then

Tð3=2; 0Þ ’ T ð1Þ ;  @T  ’ Pð1Þ ; @r  ðr;zÞ¼ð3=2;0Þ

ð1Þ

where T can be easily worked out from the given boundary conditions (since T is specified on r ¼ z þ 3=2) and P ð1Þ is an unknown to be determined. Numerical values are obtained for T by using N = 141 and N = 281 with d ¼ 0:0001.The numerical results are compared graphically with the exact solution as shown in Figs. 8 and 9. Fig. 8 shows the temperature along r = 1 ( 12 < z < 2) while Fig. 9 captures the variation of the temperature at z = 1 and z = 2 (0 < r < 2). On the whole, the numerical and exact temperature values agree well with each other.

Problem 5. Here we consider a thermal management system modeled by two homogeneous cylindrical solids as sketched in Fig. 10. The regions X1 and X2 model the computer chip and the heat sink respectively while the interface C0 (line z ¼ 0; 0 < r < r2 ) represents a thin layer of carbon nanotubes or nanocylinders of high thermal conductivity. We model the interface C0 as high conducting. A constant heat flux q0 flows into the system through z ¼ z1 ; 0 < r < r 1 . There is a uniform convective cooling at the end z ¼ z2 ; 0 < r < r2 . Elsewhere, the exterior boundary of the thermal system is thermally insulated. More precisely, the boundary conditions on the sides that are not thermally insulated are as follows:

 j1 Pðr; z1 ; 0; 1Þ ¼ q0

for 0 < r < r 1 ;

 j2 Pðr; z2 ; 0; 1Þ ¼ h½Tðr; z2 Þ  T a  for 0 < r < r 2 ; where q0 is the magnitude of the specified heat flux, h is the heat convection coefficient and T a is the ambient temperature of the system. We study the effect of the interfacial parameter a (assumed to be constant) on the thermal performance of the heat dissipation system. For this purpose, we take the radii r 1 and r2 and the lengths z1 and z2 to be r 2 =r 1 ¼ 5, r1 =z1 ¼ 1 and z2 =z1 ¼ 5.

1960

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

Fig. 8. Plots of the numerical and exact boundary temperature Tð1; zÞ for  12 < z < 2 (Problem 4).

Fig. 9. Plots of the numerical and exact boundary temperature Tðr; zÞ for z ¼ 1 and z ¼ 2 (0 < r < 2) (Problem 4).

The numerical results are obtained by employing a total of 160 elements on the exterior boundary of X1 [ X2 . To capture the thermal behaviors more accurately near the region of heating, it may be necessary to employ more elements on the side z ¼ z1 , especially near the point ðr; zÞ ¼ ðr 1 ; z1 Þ where the boundary heat flux is discontinuous. Using

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

1961

Fig. 10. A geometrical sketch of Problem 5 on the rz plane.

Fig. 11. Plots of

j1 ðT  T a Þ=ðq0 z1 Þ against z=z1 for a few selected values of a=ðj1 z1 Þ (for high conducting interfaces).

hz1 =j1 ¼ 2:5  103 and j2 =j1 ¼ 2:20, the non-dimensionalized temperature j1 ðT  T a Þ=ðq0 z1 Þ along the z-axis are plotted against z=z1 for selected values of the non-dimensionalized parameter a=ðj1 z1 Þ. In Fig. 11, the dashed line (a=ðj1 z1 Þ ¼ 0:05) approximates the plot of the non-dimensionalized temperature profile for the case in which the interface between the chip and heat sink is nearly perfectly bonded (for perfectly bonded or ideal interface,

1962

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

Fig. 12. Plots of

j1 ðT  T a Þ=ðq0 z1 Þ against z=z1 for a few selected values of a=ðj1 z1 Þ (for low conducting interfaces).

a=ðj1 z1 Þ ¼ 0). As anticipated, at a given point on the z-axis, the non-dimensionalized temperature in both the computer chip and heat sink decreases as a=ðj1 z1 Þ increases. Hence, the thin layer of carbon nanotubes or nanocylinders of high thermal conductivity enhances the heat dissipation performance of the system. Still with r 2 =r 1 ¼ 5; r 1 =z1 ¼ 1; z2 =z1 ¼ 5; hz1 =j1 ¼ 2:5  103 and j2 =j1 ¼ 2:20, we will now investigate the case whereby the interface between the chip and the sink is filled with microscopic voids. We regard this interface as low conducting. Again, we plot j1 ðT  T a Þ=ðq0 z1 Þ against z=z1 for selected values of the non-dimensionalized parameter kz1 =j1 as shown in Fig. 12. The dashed line (kz1 =j1 ¼ 100 ) gives the temperature profile for the case of a nearly ideal interface as there is negligible temperature jump across the interface at z ¼ 0. Note that the dashed lines in both Figs. 11 and 12 give the temperature profile for a nearly perfect interface. The temperature in the chip in Fig. 12 is higher than that in Fig. 11. Thus, the effect of the low conducting interface on heat flow is opposite to that of the high conducting one, that is, the low conducting interface obstructs rather than enhance the heat flow from the chip into the sink. As expected, when the obstruction of the heat flow is higher (that is, when kz1 =j1 has a lower value), the temperature jump across the interface at z=z1 ¼ 0 is bigger. Also, the differences between the temperature distributions for the different values of kz1 =j1 are much smaller in the sink compared to those in the chip. On the whole, Figs. 11 and 12 summarize the effects of the three types of interfaces – low conducting, perfectly conducting and high conducting ones – on the thermal performance of the heat dissipation system in Fig. 10. 5. Summary Boundary element procedures based on special Green’s functions are proposed for analyzing axisymmetric heat conduction across low conducting and high conducting interfaces between two dissimilar materials. As the Green’s functions satisfy the relevant interfacial conditions, the boundary element procedures do not require the interfaces to be discretized into elements, giving rise to smaller systems of linear algebraic equations to be solved. The procedures are applied to solve particular problems with known exact solutions. The numerical solutions obtained confirm the validity of the Green’s functions and the proposed Green’s function boundary element procedures. Appendix A With reference to a Cartesian coordinate system Oxyz, consider the half-spaces z < 0 and z > 0 being occupied by two dissimilar homogeneous materials having thermal conductivities j1 and j2 , respectively. We derive here Green’s functions for three-dimensional steady-state heat conduction across low and high conducting interfaces at z = 0.

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

1963

The required Green’s functions denoted by Uðx; y; z; n; g; fÞ are solutions of the partial differential equation

@2U @2U @2U þ 2 þ 2 ¼ dðx  nÞdðy  gÞdðz  fÞ; @x2 @y @z

ðA1Þ

where d denotes the Dirac-delta function. Eq. (A1) admits solutions of the form

Uðx; y; z; n; g; fÞ ¼ 

HðzÞHðfÞ þ HðzÞHðfÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ U ðx; y; z; n; g; fÞ; 4p ðx  nÞ2 þ ðy  gÞ2 þ ðz  fÞ2

ðA2Þ

where U ðx; y; z; n; g; fÞ is any solution of the partial differential equation

@ 2 U @ 2 U @ 2 U þ þ ¼ 0: @x2 @y2 @z2

ðA3Þ

Note that HðxÞ denotes the unit-step Heaviside function. The function U ðx; y; z; n; g; fÞ are given in [6] for low and high conducting interfaces for f ¼ 0. The analysis in [6] can be modified to include the case where f – 0. We take U ðx; y; z; n; g; fÞ to be of the form

8 2 3 > Z 1 < a expða uÞdu 0 3 6 7 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi5 U ðx; y; z; n; g; fÞ ¼ HðfÞ HðzÞ4 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ a1 > 2 2 2 2 2 2 0 : 4p ðx  nÞ þ ðy  gÞ þ ðz þ fÞ ðx  nÞ þ ðy  gÞ þ ðz þ f  uÞ 9 > Z 1 = expða3 uÞdu qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þHðzÞa2 > 0 ðx  nÞ2 þ ðy  gÞ2 þ ðz  f þ uÞ2 ; 8 2 3 > Z 1 < b expðb uÞdu 0 3 6 7 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi5 þ HðfÞ HðzÞ4 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ b1 > 0 : 4p ðx  nÞ2 þ ðy  gÞ2 þ ðz þ fÞ2 ðx  nÞ2 þ ðy  gÞ2 þ ðz þ f þ uÞ2 9 > Z 1 = expðb3 uÞdu qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; ðA4Þ þHðzÞb2 > 0 ðx  nÞ2 þ ðy  gÞ2 þ ðz  f  uÞ2 ; where a0 ; a1 ; a2 ; a3 ; b0 ; b1 ; b2 and b3 are constants. We assume a priori that a3 and b3 are positive constants so that the improper integrals over ½0; 1Þ in (A4) exist. It may be easily verified that (A4) is a solution of (A3) at all points ðx; y; zÞ in space. The constants a0 ; a1 ; a2 ; a3 ; b0 ; b1 ; b2 and b3 are chosen to satisfy the conditions on the interfacial conditions. A.1. Low conducting interface For the case in which z = 0 is low conducting, Uðx; y; z; n; g; fÞ is required to satisfy the interfacial conditions

    @ @ ¼ j2 ½Uðx; y; z; n; g; fÞ ; ½Uðx; y; z; n; g; fÞ @z @z z¼0 z¼0þ   @ : k½Uðx; y; 0þ ; n; g; fÞ  Uðx; y; 0 ; n; g; fÞ ¼ j2 ½Uðx; y; z; n; g; fÞ @z z¼0þ

j1

ðA5Þ

If we take a0 ¼ b0 ¼ 1, the condition on the first line of (A5) is satisfied if

 j1 a1 ¼ j2 a2 ;  j1 b2 ¼ j2 b1 :

ðA6Þ

For f < 0, the condition on the second line of (A5) is satisfied if

kða2  a1 Þ

Z

1

0

¼ j2 a2

Z 0

expða3 uÞdu k qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 ðx  nÞ þ ðy  gÞ þ ðu  fÞ 2p ðx  nÞ þ ðy  gÞ2 þ f2

1

ðu  fÞ expða3 uÞdu ½ðx  nÞ2 þ ðy  gÞ2 þ ðu  fÞ2 3=2

Using the integration by parts, we obtain

:

ðA7Þ

1964

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

Z

ðu  fÞ expða3 uÞdu 2

2

2 3=2

½ðx  nÞ þ ðy  gÞ þ ðu  fÞ 

¼

expða3 uÞ

 a3 ½ðx  nÞ2 þ ðy  gÞ2 þ ðu  fÞ2 1=2

Z

expða3 uÞdu qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : ðA8Þ ðx  nÞ2 þ ðy  gÞ2 þ ðu  fÞ2

From (A7) and (A8), it follows that

 2pj2 a2 ¼ k; j2 a2 b3 ¼ kða2  a1 Þ:

ðA9Þ

Similarly, for f > 0, we obtain

2pj2 b1 ¼ k;  j2 b1 b3 ¼ kðb1 þ b2 Þ:

ðA10Þ

Solving (A6), (A9) and (A10) gives

k ; 2pj1 k ; a2 ¼ b1 ¼  2pj2 k j2 ð1 þ Þ: a3 ¼ b3 ¼

a1 ¼ b2 ¼

j2

ðA11Þ

j1

Note that a3 and b3 are positive (as assumed). Thus, the required three-dimensional Green’s function for the case in which the interface z ¼ 0 is low conducting is given by (A2), (A4) and (A11). A.2. High conducting interface For the case in which z = 0 is high conducting, Uðx; y; z; n; g; fÞ is required to satisfy the interfacial conditions

Uðx; y; 0þ ; n; g; fÞ ¼ Uðx; y; 0 ; n; g; fÞ;

      @ @ @2    j2 ½Uðx; y; z; n; g; fÞ  j1 ½Uðx; y; z; n; g; fÞ ¼ a 2 ½Uðx; y; z; n; g; fÞ  @z @z @z z¼0þ z¼0

ðA12Þ

:

z¼0

Taking a0 ¼ b0 ¼ 1, the condition on the first line of (A12) is satisfied if

a1 ¼ a2 ; b1 ¼ b2 :

ðA13Þ

For f < 0, the condition on the second line of (A12) is satisfied if

ðj1 a1 þ j2 a2 Þ

Z

Z 0

1

ðu  fÞ expða3 uÞdu ½ðx  nÞ2 þ ðy  gÞ2 þ ðu  fÞ2 3=2 2

j1 f

þ

2p½ðx  nÞ2 þ ðy  gÞ2 þ f2 3=2 3   @2 6 1  7 expða3 uÞ 2 4qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi5 du:  @z ðx  nÞ2 þ ðy  gÞ2 þ ðu  fÞ2  0

¼ aa2

1

ðA14Þ

z¼0

Using integration by parts and the relation

3   @ 6 1 7 q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  4 5 2  @z ðx  nÞ2 þ ðy  gÞ2 þ ðz  f þ uÞ2  2

2

2

z¼0

3

2

¼

@ 6 1 7 4qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi5; @u2 2 2 2 ðx  nÞ þ ðy  gÞ þ ðu  fÞ

ðA15Þ

we obtain

Z 1 ðj1 a1 þ j2 a2 Þ expða3 uÞdu q qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi þ a  ð j a þ j a Þ 3 1 1 2 2 2 2 2 3=2 2 2 2p½ðx  nÞ þ ðy  gÞ þ f  0 ðx  nÞ þ ðy  gÞ þ f2 ðx  nÞ2 þ ðy  gÞ2 þ ðu  fÞ2 8 9 > > Z 1 < = f a3 expða3 uÞdu 2 q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ a  ¼ aa2  : ðA16Þ 3 2 2 3=2 2 > > 0 : ½ðx  nÞ þ ðy  gÞ þ f  ðx  nÞ2 þ ðy  gÞ2 þ f2 ðx  nÞ2 þ ðy  gÞ2 þ ðu  fÞ2 ;

j1 f

From (A16), it follows that

E.L. Chen, W.T. Ang / Applied Mathematical Modelling 37 (2013) 1948–1965

j1 ¼ 2paa2 ; ðj1 a1 þ j2 a2 Þ ¼ aa2 a3 :

1965

ðA17Þ

Similarly, for f > 0, we obtain

j2 ¼ 2pab2 ; ðj1 b2 þ j2 b1 Þ ¼ ab2 b3 :

ðA18Þ

Solving (A13), (A17) and (A18) gives

a1 ¼ a2 ¼ 

j1 2pa

;

j2 ; b1 ¼ b2 ¼  2pa 1 a3 ¼ b3 ¼ ðj1 þ j2 Þ: a

ðA19Þ

Note that a3 and b3 are positive (as assumed). Thus, the required three-dimensional Green’s function for the case in which the interface z = 0 is high conducting is given by (A2), (A4) and (A19). References [1] Y. Benveniste, On the decay of end effects in conduction phenomena: a sandwich strip with imperfect interfaces of low or high conductivity, J. Appl. Phys. 86 (1999) 1273–1279. [2] A. Desai, J. Geer, B. Sammakia, Models of steady heat conduction in multiple cylindrical domains, J. Electron. Packag. – Trans. ASME 128 (2006) 10–17. [3] J.R. Berger, A. Karageorghis, The method of fundamental solutions for heat conduction in layered materials, Int. J. Numer. Methods Eng. 45 (1999) 1681– 1694. [4] W.T. Ang, K.K. Choo, H. Fan, A Green’s function for steady-state two-dimensional isotropic heat conduction across a homogeneously imperfect interface, Commun. Numer. Methods Eng. 20 (2004) 391–399. [5] W.T. Ang, Non-steady state heat conduction across an imperfect interface: a dual-reciprocity boundary element approach, Eng. Anal. Boundary Elem. 30 (2006) 781–789. [6] X. Wang, L.J. Sudak, 3D Green’s functions for a steady point heat source interacting with a homogeneous imperfect interface, J. Mech. Mater. Struct. 1 (2006) 1269–1280. [7] C.A. Brebbia, J.C.F. Telles, L.C. Wrobel, Boundary Element Techniques, Theory and Applications in Engineering, Springer-Verlag, Berlin/Heidelberg, 1984.