Heat conduction in a semi-infinite medium with a spherical inhomogeneity and time-periodic boundary temperature

International Journal of Heat and Mass Transfer 55 (2012) 618–628

Contents lists available at SciVerse ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Heat conduction in a semi-infinite medium with a spherical inhomogeneity and time-periodic boundary temperature A. Rabinovich ⇑, G. Dagan, T. Miloh School of Mechanical Engineering, Tel-Aviv University, Tel-Aviv 69978, Israel

a r t i c l e

i n f o

Article history: Received 17 May 2011 Received in revised form 5 October 2011 Accepted 13 October 2011 Available online 24 November 2011 Keywords: Heat conduction Time-periodic Semi-infinite medium Heterogeneous medium Analytical solution Perturbation expansion

a b s t r a c t We solve the problem of heat conduction in a homogeneous media below a planar boundary subjected to time-periodic temperature (of frequency x), in the presence of a spherical inhomogeneity (of radius R), whose center is at distance d > R from the boundary. In the absence of the sphere, the well known one dimensional pffiffiffiffiffiffiffiffiffiffiffiffiffi solution can be regarded as an oscillating thermal boundary layer of displacement thickness d ¼ 2a=x, where a is the heat diffusivity. The general solution depends on four dimensionless parameters: d/R, d/R, the heat conductivity ratio j and the heat capacity ratio C. An analytical solution is derived as an infinite series of Bessel functions, which converges quickly. The results are illustrated and analyzed for a given accuracy and for a few values of the governing parameters. The general solution can be simplified considerably for asymptotic values of the parameters. A first approximation, obtained for R/d  1, pertains to an unbounded domain. A further approximate solution, for R/d  1, while j and C are fixed, can be regarded as pertaining to a quasi-steady regime, and is similar in structure to Maxwell’s solution for steady state. However, its accuracy deteriorates for j  1, and a solution, coined as the insulated sphere approximation, is derived for this case. Comparison with the exact solution shows that these approximations are accurate for a wide range of parameter values. Besides providing insight, they can be employed for solving in a simple manner more complex problems, e.g. effective properties of a heterogeneous medium made of an ensemble of spherical inclusions. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction We consider here the problem of heat conduction into the half space, below a planar boundary subjected to time-periodic temperature. The medium of constant properties (conductivity, diffusivity) contains a spherical inclusion of different constant properties. While the solutions for a homogeneous medium or an isolated sphere subjected to given time dependent temperature on its boundary are well documented ([1] – Chapters 7, 8 and 11), much less is known for configurations similar to the present one. There is vast literature on the steady state problem (applicable to similar processes like electrical conduction, diffusion and flow through porous media) subjected to a uniform field at infinity, starting with the classical Maxwell solution [2]. Among the few previous works on the unsteady problem, the recent studies [3,4] are most relevant to the one pursued here. Though we follow their general approach of expanding the solution of the governing Helmholtz equation in a series of eigenfunctions, we explicitly derive the solution for the half space with periodic temperature variations, whereas they provide a general approach for an infinite domain solely. Furthermore, ⇑ Corresponding author. Tel.: +972 3 640 8227; fax: +972 3 640 7334. E-mail address: [email protected] (A. Rabinovich). 0017-9310/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2011.10.049

we present here for the first time a full analytical solution to the above problem together with simple approximate solutions which lend themselves to physical interpretation. Such solutions may be extremely useful when solving more complex problems. The recent article of [5] also addresses a similar time dependent problem for a spherical inclusion, yet it is limited to a highly conducting sphere with only radial dependence of temperature and it assumes heat generation. Another article which solves a mathematical problem similar to that of this work and presents some approximations as well is [6]. Nevertheless, the problem considered there is of diffusive interaction between two ideal sinks, which is not time periodic and has simpler boundary conditions. The applications of primary interest to us are of a geophysical nature, e.g. the impact of heterogeneity upon heat flow through the earth crust under periodic diurnal or annual temperature variations or the related problem of fluid flow through elastic porous formations in the soil. However, other potential applications like heat exchange between blood tissue and embedded blood vessels [7], the hydrocooling of fruits or vegetables of spherical shape [8] and experimental methods for specifying the thermal diffusivity of materials [9], can also be envisaged. Similarly, the present solution may constitute the first step toward determining equivalent properties of heterogeneous media under unsteady conditions.

619

A. Rabinovich et al. / International Journal of Heat and Mass Transfer 55 (2012) 618–628

Thus, we consider the problem addressed here as of fundamental nature in view of the possible ramifications of the solution. The plan of the paper is as follows: the problem is stated mathematically in Section 2 and its exact solution is derived subsequently in Section 3, expressed in terms of the various dimensionless parameters of the problem. Two simple solutions coined here as complex Maxwell (CM) and the insulated sphere approximations (ISA) are derived in Section 4 and their accuracy is determined by comparison against the exact one. The impact of the spherical inhomogeneity upon the surface heat flux, which is perpendicular to the planar boundary, is examined in Section 5. The paper is concluded with a summary and discussion in Section 6.

We consider a semi-infinite domain subjected to periodic temperature variation and enclosing a spherical body of radius R (Fig. 1). The boundary temperature is assumed to be constant in space and to vary sinusoidally in time. We further assume a medium of constant ambient thermal diffusivity, aex, with a spherical inhomogeneity of radius R and diffusivity  ain located  at a depth , where K is d P R from the surface (aex=in ¼ K ex=in = qex=in  C ex=in p the heat conductivity, q the density and Cp the specific heat of the two media). For convenience we use both cartesian and spherical coordinates defined by x = r sin h cos u, y = r sin h sin u and z = r cos h, taking the origin of the axis to be at the sphere center (see Fig. 1). The governing equations for the temperature of the sphere interior (denoted ‘‘in’’) and the exterior temperature (denoted ‘‘ex’’) are

z6d

ð1Þ

ð2aÞ ð2bÞ

where x is the prescribed frequency of the surface temperature and Thom is the solution for the temperature field in the homogeneous medium of diffusivity aex. The usual conditions for continuity of temperature and heat flux across the spherical boundary are as follows:

T ex ¼ T in ; r ¼ R @T ex @T in ¼j ; r¼R @r @r

ð4Þ

Tb ex ðx; y; dÞ ¼ 1;

Tb ex ! Tb hom

ðfor r ! 1Þ

ð5Þ

The solution of the heat flow problem in a homogeneous media was presented and analyzed in [1] (Section 2.6) and is given by

pffiffiffiffiffiffiffi Tb hom ¼ ekex ðzdÞ ; i:e: j Tb hom j ¼ eðzdÞ= 2aex ¼ eðzdÞ=d pffiffiffiffiffiffiffiffiffi argð Tb hom Þ ¼ ðz  dÞ= 2aex ¼ ðz  dÞ=d

ð6Þ

Rd ffi b hom jdz ¼ pffiffiffiffiffiffiffiffiffiffi where d ¼ 1 j T 2aext is the displacement thickness of b hom j, a convenient length the thermal boundary layer defined by j T scale characterizing the depth of penetration into the medium. It is convenient to recast the problem for T 0ex , the normalized perturbation temperature associated with the spherical inclusion, defined by

  Tb ex ¼ ekex d T 0ex þ Tb hom ¼ ekex d T 0ex þ ekex z ;

rP1

ð7Þ

while T 0in is defined by

with the boundary conditions

T ex ¼ T 0 cosðxtÞ; z ¼ d T ex ¼ T hom ; r ! 1

r2 Tb ex=in ¼ k2ex=in  Tb ex=in

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi where kex=in ¼ ði þ 1Þ= 2aex=in . The boundary conditions (2a,b) satisb ex are consequently fied by T

2. Mathematical statement

@T ex=in ¼ aex=in  r2 T ex=in @t

where j = Kin/Kex. We reformulate the problem with the aid of dimensionless variables, to simplify its parameter dependence. Thus, variables are normalized by using R as length, x1 as time and T0 as temperature scales. Hence, the dimensionless variables are defined by x/R, y/R, z/R, r/R, d/R, T/T0, aex/in/xR2, xt and for convenience we maintain the original notation for these variables in the following. The boundary condition (2a) implies the solution is of the form b eit  allowing us to solve for the time independent variable, T ¼ Re½ T b . The following modified Helmholtz equations are then derived T b ex=in from (1) for the complex temperatures T

Tb in ¼ ekex d T 0in ;

r61

ð8Þ

Thus, it is found from (4), (5) and (3a,b) that the complete set of equations and boundary conditions satisfied by T 0ex=in are as follows:

r2 T 0ex  k2ex  T 0ex ¼ 0;

z 6 d; r P 1

ð9aÞ

r2 T 0in  k2in  T 0in ¼ 0;

z 6 d; r 6 1

ð9bÞ

T 0ex ¼ 0;

z¼d

ð10aÞ

ð3aÞ

T 0ex ¼ 0;

r!1

ð10bÞ

ð3bÞ

T 0in  T 0ex ¼ ekex z ;

j

@T 0in @r



@T 0ex @r

¼

r¼1 @ kex z e ; @r

ð11aÞ r¼1

ð11bÞ

It is seen that for j = 1 and aex = ain the solution of the system bT b hom : (9)–(11) is given by T 0ex  0; T 0in ¼ ekexz , i.e. T  in Defining the ratio C ¼ qex C ex p = qin C p

along with the rela-

tionships aex = d2/2 and ain = jCd2/2, renders the temperature fields as functions of four independent dimensionless parameters: j, d, C and d. 3. Exact solution The general solution of (9a,b) for T0 is given for the sphere exterior and interior, respectively, by [10]: 1 X n X

Anm eimu Pm n ðlÞ

kn ðkex  rÞ ; kn ðkex Þ

rP1

ð12aÞ

Anm eimu Pm n ðlÞ

g n ðkin  rÞ ; g n ðkin Þ

r61

ð12bÞ

n¼0 m¼n

Fig. 1. The setup of the problem as formulated by (1) and (2a,b). The sphere and its image are denoted by a solid and dashed circle respectively and a planar boundary is located at a distance d above the sphere center.

1 X n X n¼0 m¼n

620

A. Rabinovich et al. / International Journal of Heat and Mass Transfer 55 (2012) 618–628

where gn and kn are the regular and singular modified spherical Bessel functions, respectively, related to the modified Bessel functions pffiffiffiffiffi p K 1 ðkrÞ and of order equal to half an odd integer: kn ðkrÞ ¼ 2kr nþ2 pffiffiffiffiffi p g n ðkrÞ ¼ 2krInþ1 ðkrÞ. Anm and Bnm are complex coefficients, Pm n are

where the factor (n + 1/2) appears for convenience only and it is repffiffiffiffiffiffiffi minded that kex = (i + 1)/d and kin ¼ ði þ 1Þ=ðd jC Þ. Simple geometry bears the relation between the two coordinates:

Legendre polynomials and l = cos h. To satisfy the surface boundary condition, (10a), we employ the image technique which implies supplementing T 0ex=in , defined in the entire space with boundary

ð15Þ

2

condition T 0ex ! 0 for r ? 1 (denoted T 0sph ex ), with the temperature 0sph 0 0 field T 0im ex ¼ T ex ðr ; h Þ associated with a fictitious sphere located symmetrically at a distance d above the surface (Fig. 1). Thus 0im 0im T 0ex ¼ T 0sph ex þ T ex where T ex acts as an exterior forcing term and Eqs. (11a,b) have to be replaced by 0 kex z 0 T 0in  T 0sph  T 0sph ex ¼ e ex ðr ; h Þ;

r¼1

@T 0 @T 0sph @ @T 0sph ðr0 ; h0 Þ ; j in  ex ¼ ekex z  ex @r @r @r @r

ð13Þ r¼1

the coordinate system (r0 , h0 ) having the origin at z = 2d, situated at the center of the image (Fig. 1). This results in the fulfillment of the condition T 0ex ¼ 0 on z = d, b hom ¼ 1 there. By using leaving only the homogeneous solution T this technique and due to symmetry in the u coordinate, we arrive at the following general solution of (9a,b) with conditions (10a,b) 1 X kn ðkex  rÞ ðn þ 1=2ÞAn Pn ðlÞ kn ðkex Þ n¼0 1 X kn ðkex  r 0 Þ ðn þ 1=2ÞAn  Pn ðl0 Þ; r P 1 kn ðkex Þ n¼0 1 X g ðkin  rÞ ðn þ 1=2ÞBn n T 0in ¼ Pn ðlÞ; r 6 1 g n ðkin Þ n¼0

T 0ex ¼

ð14aÞ ð14bÞ

r0 ¼

3.1. Results The solution T0 (r, h) is computed at any point in space by truncation of the series in (14a,b). The number of terms 0 6 N < 1 used for computation is determined by the desired accuracy in satisfying the boundary conditions (11a,b). The results displayed in the  h i b b  b following are for a 0.5% accuracy, i.e. j T < 0:005 ex jj T in j=j T ex j r¼1  h i  b  b in =@rj=j@ T b ex =@rj < 0:005. Some examples and j@ T ex =@rj  jj@ T r¼1

for values of N necessary to reach this accuracy are given in Table 2. For illustration, Fig. 2 shows the isotherms for a cross section of four spheres with different conductivity ratios (j) and depths (d). b , related to The plots are for t = 0 and they represent Tjt¼0 ¼ Re½ T T0 through (7) and (8). While j has a direct impact on T 0in via kin in (14b), d impacts T 0ex=in only indirectly due to the presence of the image. It is seen that the behavior of the temperature field inside the sphere is different depending on whether j > 1 or j < 1. Thus, for j > 1 the temperature variation inside the sphere is relatively small and as j ? 1 the temperature becomes constant. When j < 1 the temperature variation within the sphere is enhanced and for j ? 0 (insulated sphere) the temperature changes

(c) 3

2

2

1

1

0

0

z

3

−1

−1

−2 −3 −3

−2

−1

(b)

0

x

1

2

3

0.8

−2

0.6

−3 −3

0.4

−2

−1

−2

−1

(d)

x

0

1

2

3

0

1

2

3

1.5

1

0.2

1

0

0

0

z

z

1.5

−1 −2 −3 −3

l0 ¼ ð2d  rlÞ=r0

The coefficients An and Bn are subsequently determined by applying (11a,b) and solving a linear set of equations (the full derivation is given in Appendix A).

(a)

z

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 r 2 þ 4d  4rdl;

−1 −2

−2

−1

0

x

1

2

3

−3 −3

x

b , of a cross section through the center of spheres of different conductivities and depths. (a) and (b) depict spheres that have a larger Fig. 2. The temperature, Tjt¼0 ¼ Re½ T conductivity than their surroundings, j = 5, at depths d = 3 (a) and d = 1.5 (b) from the surface. (c) and (d) correspond to j = 0.2 for d = 3 (c) and d = 1.5 (d). d = 2 and C = 1 in all plots.

621

A. Rabinovich et al. / International Journal of Heat and Mass Transfer 55 (2012) 618–628

(a)

(b)

b j, as a function of z for the vertical line passing through the center of the sphere. Parameter values are d = 2 (a) and d = 10 Fig. 3. The amplitude of the complex temperature, j T (b) while C varies as shown in the figure legends. j = 1 and d = 3 for both figures. Notice that the C = 100, C = 1 and C = 0.1 lines in (b) coincide.

drastically tending to zero. As for the impact of d, it is seen that the isotherms are horizontal near the surface boundary for the large d = 3, whereas the curved isotherms of the perturbation reach the boundary for d = 1.5 when the condition (10a) affects the solution. b eixt  ¼ It is emphasized that in the general expression T ¼ Re½ T b j cosðxt þ cÞ where c ¼ argð T b Þ, both j T b j and c depend on x, y, z jT and the parameters of the problem. In the following we focus on b j, while the phase c can be determined by a similar illustrating j T procedure. b , while j = 1, Fig. 3 deTo illustrate the impact of C and d on T b j as a function of z for the vertical line passpicts the amplitude j T ing through the center of the sphere. A large d (Fig. 3b) results in a b hom j corresponding to the C = 1 line in the plot, practically linear j T b hom j drops exponentially from z = d. while for a small d (Fig. 3a), j T The change brought on by the parameter C is shown for j = 1 so b hom . It is reminded that the steady it is most easily compared to T state problem has no dependence on the heat capacity of the two media and for j = 1 the inhomogeneity is no more present. Here, raising C above unity is of meager affect while values below unity have a much more noticeable impact. When C ? 0 it is seen in b ðr 6 1Þ, tends toward zero, as displayed in the (14b) that Tin, i.e. T plotted amplitude profiles.

4. Approximate solutions The complete solution of T0 depends on four independent parameters and is quite involved. We present and analyze now a few approximations of T0 for particular ranges of parameter values, toward obtaining simple algebraic expressions for later use in more complex problems. Moreover, these approximations serve as a useful tool in grasping the physics of the problem. Finally, we show in Appendix C, that some of the approximations can be obtained by a perturbation expansion of the general equations, thus circumventing the need to solve exactly first and enabling the use of similar approximations in more complex cases where an analytical solution is not available.

4.1. A large depth approximation If the sphere is far enough from the boundary, i.e. d  1, the effect of the image sphere, i.e. the second term in (14a), can be neglected. This is equivalent to solving for T0 for a sphere in an infinite domain so that the boundary condition (10a) at z = d is immaterial, while the ambient temperature satisfies the boundary b hom ¼ 1 at z = d. The coefficients of the expansion of condition T T 0ex=in , pertaining to (14a,b), are now given by the closed form expressions

 R 1 pffiffiffiffiffiffiffiffiffi j=C  g~m  l ekex l Pm ðlÞdl 1 Am ¼ pffiffiffiffiffiffiffiffiffi ~m  j=C  g~m k Z 1 ekex l P m ðlÞdl Bm ¼ Am þ

ð16aÞ ð16bÞ

1

~m ¼ 1 where k km ðkex Þ

h

@km ðxÞ @x

i x¼kex

1 and g~m ¼ gm ðk in Þ

h

@g m ðxÞ @x

i x¼kin

(see Appendix

B). While the coefficients in (16a,b) can be easily calculated numerically, they have analytical expressions which are obtained by replacing the integrals with the expressions given in Appendix A (see (A8)). Since T 0ex=in (14a,b) is independent of d, the temperature b is influenced by d in a simple manner through the complex conT stant ekex d in (7) and (8). Next, we compare the exact solution with this approximation, given by (16a,b), in order to determine the range of depths for which it is accurate. Adopting for the purpose of comparison the most stringent condition, we select the point (r = d, h = 0) on the upper boundary, which is closest to the image. We find that for d = 3 the approximate solution differs from the exact one at this point by less than 5%. This is true for any choice of parameters, except for extremely small values of C not considered here. As a matter of fact, for most parameter values the difference is less than 1%. In the rest of this Section we assume d > 3 and refer to the large depth solution solely. An additional limitation regarding d and related to the depth d stems from the requirement that the sphere is not outside the

622

A. Rabinovich et al. / International Journal of Heat and Mass Transfer 55 (2012) 618–628

b hom  1, where the temperboundary layer, in a region in which T ature is practically constant and negligible. Thus, we impose the b hom j at the center of the sphere following condition for j T

ed=d >  ! d > 

4.3. The insulated sphere approximation

d lnðÞ

ð17Þ

For example, selecting  = 0.05 and d = 3 yields d > 1. A smaller value of d will be regarded as leading to a degenerate solution. Though the solution for d > 3 (16a,b) is less complicated than the exact one, even simpler approximations, which lend themselves to physical interpretations, are derived in the sequel by further limiting the range of values of the other parameters. 4.2. The complex Maxwell (quasi-steady) approximation We consider next the limit d pffiffiffiffiffiffi ffi 1, or jkexj  1, while C and j are given, such that jkin j ¼ jkex j= C j  1 as well. From a physical perspective this can be regarded as the case of a medium of given properties, at the low frequency limit x  1. Indeed, the dimensional diffusivities aex/in were normalized with respect to x, i.e. jkex/inj with respect to x1/2. Thus, we can regard this approximation as a quasi-steady one. In the same vein, the limit d  1 implies that b hom can be approximated by a linear profile in the boundary conT ditions (11a,b) obeyed by T0 , i.e. the inclusion finds itself in a field of constant temperature gradient. To obtain the approximate solution pertinent to this limit, the solution based on (16a,b) is expanded for jkex/inj ? 0 and first-order terms are retained. With details given in Appendix B, the final result is

T 0M;ex ¼ kex

 1j l ; 2 þ j r2

respectively. An approximation pertinent to this configuration is examined next.

T 0M;in ¼ 1 þ kex

3 lr 2þj

ð18Þ

It is seen that this approximate solution resembles the classical Maxwell solution for the temperature field in steady state [2], obtained for TM which satisfies the Laplace equation, with the boundary condition TM ? 1 + kexz = 1 + kexlr for r ? 1. The direct derivation of (18) by a perturbation expansion of the exact Eqs. (9a,b), (10b) and (11a,b) is discussed in Appendix C. It is emphasized that unlike the Maxwell solution, T 0M is complex, and this makes an important difference when substituting (18) in (7) and b ex=in . Indeed this leads to (8) to obtain T

  1j z kex z Tb M;ex ¼ ekex d ðT 0M;ex þ ekex z Þ ¼ ekex d kex ð19aÞ þ e 2 þ j r3

  1j z þ ekex z Tb M;in ¼ ekex d ðT 0M;in  1  kex z þ ekex z Þ ¼ ekex d kex 2þj ð19bÞ as the temperature field, with kex = (i + 1)/d, which still has to be multiplied by eit in order to arrive at the unsteady TM (1–3). b M (19a,b) is an inconsistent approximaIt is emphasized that T b hom ¼ ekex ðzdÞ in (7) at first order in tion as we did not expand T kex. Furthermore, the interior temperature in (18) incorporates the linear background, 1 + kexz, and therefore we replaced this value with ekex z in (19b) to maintain consistency between the interior b M with and exterior. This has a beneficial effect when comparing T the exact solution, as shown in the sequel. At any rate, even if we b M eit is difreplace ekex z by 1 + kexz in (19a,b), the real part of T M ¼ T ferent from the Maxwell steady state solution due to the difference in phase among the different terms and for this reason we have coined (18) as the complex Maxwell pffiffiffiffiffiffiffi solution. Finally, the limit jkin j ¼ jkex j= C j  1 underlying T 0Min , implies that the accuracy of the CMA (complex Maxwell approximation) deteriorates as j or C decrease, i.e. if the conductivity Kin of the sphere or the heat capacity qex  C ex p of the ambient medium are much smaller than those of the exterior and interior of the sphere,

This time we expand the large depth solution T0 of (16a,b) for ffi pffiffiffiffiffiffi jkexj ? 0 and j ? 0, while jkin j ¼ jkex j= C j ¼ 0ð1Þ and retain the leading order terms in the solution (details in Appendix B). The result is

T 0IS;ex ¼

kex l ; 2 r2

T 0IS;in ¼

g 0 ðkin rÞ 3 g ðkin rÞ l; þ kex 1 g 0 ðkin Þ 2 g 1 ðkin Þ

rP1

ð20aÞ r61

ð20bÞ

where g0(kin  r)/g0(kin) = 1/r  sinh(kinr)/sinh(kin) and g1(kin  r)/ g1(kin) = (1/r)[kin cosh(kinr)  sinh(kinr)]/[kin cosh(kin)  sinh(kin)]. It is seen that the exterior field (20a) is precisely the complex Maxwell solution (18) for j = 0, corresponding to a quasi steady state and an insulated sphere. In contrast, the interior field is the solution of the Helmholtz interior equation (9a) with a boundary condition of given temperature (3kexl)/2 induced by the exterior field (20a) on the sphere. A similar solution was obtained by [1] (Section 9.3, Eq. (12)) as a particular case of a periodic given temperature on the sphere surface. The insulated sphere approximation converges to [1] for Kex ? 1, when the exterior temperature is transmitted instantaneously to the sphere surface. From a physical perspective this approximation can be viewed as the one pertaining to steady state and linear temperature profiles at infinity at the limit of Kex  Kin, i.e. an ambient medium of much larger conductivity than that of the sphere. As a result, the exterior temperature is quasi-steady, whereas the unsteady interior field results from satisfying the temperature continuity condition (11a) solely. 4.4. Accuracy of the approximations We now analyze the validity of the above approximations by comparing them to the exact solution. The comparison is carried b j out by calculating the normalized difference, v~r ¼ jj T b approximate jj=j T b j at every point in space, then choosing jT v ¼ maxðv~r Þ to represent the relative error. Some values of v for specific cases are given in Table 2. v should be chosen according to the desired accuracy, hence defining a range of parameter values for which the approximations are adequate. We have previously discussed that d P 3 is a criterion which yields v < 5% for the large depth approximation. In Fig. 4, we test the validity of the two other approximations when d = 3. This is carried out in the plane (d, j) which are regarded as independent variables. Each contour represents the v = 5% line for a certain approximation, thus designating regions of values of j and d which render the approximation accurate. The parameter C is taken one in this figure, yet similar analysis was done for C = 0.1 and C = 10 with results summarized in Table 1. Fig. 5 depicts similar contour lines to the previous figure, but for varying d. When d < 3 the boundary condition (10a) can no longer be neglected and the temperature field associated with the sphere Table 1 Thumb rules for d P 3. Complex Maxwell approximation

Insulated sphere approximation

Degenerate solution (sphere is outside boundary layer)

C=1 C = 10

d > 3, j > 0.15 d > 3, arbitrary j

d < d/3 d < d/3

C = 0.1

d > 7, j > 0.2

d > 3, j < 0.15 d > 3, j < 0.15 (superflous) d > 7, j < 0.2

d < d/3

A. Rabinovich et al. / International Journal of Heat and Mass Transfer 55 (2012) 618–628

image impacts the result. To account in an approximate manner for this effect in the region close to the boundary z = d, we propose a   simple modification to the Maxwell T 0M;ex and insulated sphere   T 0IS;ex exterior fields, which are doublet singularities at r = 0, by supplementing them with an identical doublet at the image center r0 = 0. Thus (18) and (20a) become

T 0M;ex ¼ kex

  1j l l0 ;  2 þ j r2 r02

rP1

ð21Þ

and

T 0IS;ex ¼

1 l l kex ð 2  02 Þ; 2 r r 0

rP1

ð22Þ

623

respectively, where (r0 , l0 ) are shown in Fig. 1 and given in (15). This modification does not take into consideration the interaction between the sphere and its image and therefore as the sphere approaches the surface (d ? 1) the approximations are less accurate (see Fig. 5). Nevertheless, the approximation proves to be excellent up to a small distance, seen in Fig. 5 by the nearly identical contours of d = 3 and d = 1.5. 4.5. Discussion of results Inspection of Figs. 4 and 5 reveals that the simple Maxwell and insulated sphere approximations apply to large ranges of the problem parameters. Generally, the (d, j) plane can be divided into 3 regions, for fixed values of C and d: (i) for small d < d1(j) (which still

Fig. 4. Contours of relative error v = 0.05 both the complex Maxwell (19a,b) and insulated sphere approximations for d = 3 and C = 1. The latter  as a function  of d and j for   b IS;ex ¼ ekex d T 0 þ ekex z and T b IS;in ¼ ekex d T 0  1  kex z þ ekex z . Notice that the y-axis is logarithmic and v < 0.05 to the right of the contours. approximation is given by T IS;ex IS;in The points marked (a)–(f) correspond to the subplots in Fig. 6.

Fig. 5. Same as Fig. 4 but for varying d as shown in the figure legend. For d = 3 the approximations are the same as in the previous figure, yet for d = 1.5 and d = 1.1 the modified approximations imply replacing the exterior perturbations with Eqs. (21) and (22).

624

A. Rabinovich et al. / International Journal of Heat and Mass Transfer 55 (2012) 618–628

obey Eq. (17)), there is need to employ the exact solution. This reb hom j near the sphere or besults either from the steep variation of j T cause of the impact of the image; (ii) for d > d1 and j > j2(d) the

CMA is accurate and (iii) for d > d1 and j < j3(d) the ISA proves to be accurate. A positive finding is that the last two regions overlap, allowing for the use of the simple approximations in a large do-

(a)

(b)

(c)

(d)

(e)

(f)

b; T b homogeneous (6), the CMA (19a,b) and ISA (caption of Fig. 4) as a function of z on the vertical line passing through the sphere center. Each figure is for Fig. 6. The amplitude of T a different value of d and j corresponding to the lettered points marked in Fig. 4 and as follows: (a) d = 1, j = 10, (b) d = 5, j = 10, (c) d = 2, j = 0.01, (d) d = 6, j = 0.01, (e) d = 1.5, j = 0.1, (f) d = 2.2, j = 0.56. C = 1 and d = 3 in all figures.

^j=jq ^homogeneous j  1, on the upper boundary calculated by the exact solution, the complex Maxwell approximation and the large depth Fig. 7. The heat flux perturbation: jq approximation for C = 1, d = 10, d = 1.5 and values of j shown in the figure legend. The image effect for the approximations is taken into consideration simply by taking twice the perturbation heat flux in a similar manner to (21) and (22).

625

A. Rabinovich et al. / International Journal of Heat and Mass Transfer 55 (2012) 618–628

Table 2 Values of the temperature at the point on the sphere closest to the surface boundary (r = 1, h = 0) and the surface heat flux directly above the sphere (r = d, h = 0, also x = 0 in Fig. 7) b j=j T b hom j  1Þ, (c) is the heat flux ðjq ^j=jq ^homogeneous j  1Þ, (d) is the for different parameters: d, j, and C shown in the table and d = 2d. Details are as follows: (b) is the temperature ðj T b M j=j T b j  1Þ, (e) is the error (v) of the ISA ðj T b IS j=j T b j  1Þ and (a) is the number of terms (N) in the series solution (14a,b) used for computation for an error (v) of the CMA ðj T admissible error in temperature of 0.005 (see Section 3.1). d=3

d=2

d = 1.5

d = 1.1

d=3

d=2

d = 1.5

d = 1.1

j = 10, C = 0.1

(a) (b) (c) (d)

3 0.137 0.1545 0.0259

3 0.2233 0.5451 0.0936

5 0.3036 1.432 0.1969

12 0.3188 6.585 0.3393

j = 10, C = 2

(a) (b) (c) (d)

3 0.1131 0.0585 0.0114

4 0.1563 0.2204 0.0209

6 0.185 0.6195 0.0578

12 0.1666 2.9815 0.1784

j = 10, C = 10

(a) (b) (c) (d)

3 0.1125 0.0544 0.0114

3 0.1545 0.2061 0.02

5 0.1819 0.5824 0.0559

12 0.1630 2.8153 0.1747

j = 2, C = 0.1

(a) (b) (c) (d)

3 0.0699 0.1115 0.0319

3 0.1335 0.3727 0.097

5 0.1846 0.8851 0.1918

12 0.1294 2.6652 0.3575

j = 2, C = 2

(a) (b) (c) (d)

3 0.0381 0.0158 0.0038

3 0.0515 0.0591 0.0069

4 0.0574 0.1608 0.0199

9 0.0365 0.5772 0.064

j = 2, C = 10

(a) (b) (c) (d)

3 0.0372 0.0117 0.0038

3 0.0494 0.0451 0.0061

4 0.0542 0.1266 0.018

9 0.0341 0.4758 0.0596

j = 0.1, C = 0.1

(a) (b) (c) (d) (e)

3 0.0147 0.0406 0.2538 0.0671

3 0.0028 0.0593 0.9843 0.1671

4 0.0088 0.042 2.53 0.33

9 0.008 0.0633 7.8055 1.0524

j = 0.1, C = 2

(a) (b) (c) (d) (e)

4 0.071 0.0405 0.0073 0.011

5 0.0973 0.1423 0.0231 0.016

8 0.1027 0.3476 0.0629 0.0557

12 0.0409 0.7759 0.1933 0.2129

j = 0.1, C = 10

(a) (b) (c) (d) (e)

4 0.0721 0.0446 0.0076 0.01

6 0.0997 0.1566 0.0241 0.0166

8 0.1059 0.3817 0.0648 0.0575

14 0.0423 0.8475 0.1964 0.2160

main of the plane. Some thumb rules for using the approximations are given in Table 1. It is emphasized that the regions are dependent on the values of C and d (for d < 3) and on the choice of v. Adopting a different value of v will generate similar regions yet with parallel shifted contours. For illustration, Fig. 6 depicts temperature amplitude profiles along the vertical line passing through the sphere center. The exact b jÞ and both approximations are plotted for parameter solution ðj T values matching the corresponding points marked in Fig. 4. It is seen (point a) that for j > 1 and small d, the homogeneous solution is far from linear and therefore the CMA does not hold, though it is much closer to the exact one than the ISA, with the maximal deviation on the sphere surface. In contrast, for a large d, the point (b) belongs to the zone in which the CMA is indeed accurate. For small d and j  1 (point c) both approximations are not applicable, though the ISA is closer to the exact one, reaching the largest deviation at the bottom of the sphere while the CMA deviates inside the sphere. Increasing d for the same value j  1 as in (c), (point d) the ISA becomes accurate. The case of j < 1 but larger than (c, d) and a small d (point e) the behavior is similar to (c), while for a large d (point f) the CMA applies while an ISA fails, and a deviation at the top of the sphere is observed. 5. Heat flux on the upper boundary The effect of the presence of the sphere on the heat flux through the upper boundary is of interest in some applications. Thus, in a geophysical context, measurements of the heat flux may provide, by an inverse procedure, information about the presence of a body of different properties in the subsurface. The heat flux normal to the upper boundary is obtained by difb =@z . The heat flux pertur^ ¼ ½K ex @ T ferentiation of (7) and (8), q z¼d ^j=jq ^homogeneous j  1, is depicted in Fig. 7 for changing values bation, jq of j. A sphere which is more (less) conductive than the surrounding media will increase (decrease) the downward heat flux. Approximations of the surface heat flux are also shown in Fig. 7 with degree of accuracy similar to that shown in Fig. 4. Computations show that as j grows (diminishes) so does the heat flux

perturbation until it converges to a finite value for a perfectly conducting (non-conducting) sphere. The largest perturbation (120%) is at x = 0 and values of the heat flux perturbation for this point are given in Table 2. It is seen that as d nears unity the effect of the sphere on the heat flux is substantial. Furthermore, as shown in Table 2, a larger value of C does not have much effect on the results, yet as C tends towards zero the heat flux anomaly grows subb in j ! 0 and stantially. This is also observed in Fig. 3, where j T b j from unity at the surface consequently there is a drastic drop in j T to zero inside the sphere, thus causing an increase in the heat flux perturbation. 6. Summary and discussion We present an analytical solution to heat conduction near a spherical inhomogeneity for a half space of given harmonic temperature on its boundary. The solution is valid at any point in space and time and for all parameter values. It is obtained by using the general solution to the Helmholtz equation, an image technique to implement the boundary condition and the addition theorem for the Helmholtz equation. We derive two simple approximate solutions when the ambient temperature amplitude profile becomes close to linear near the sphere and when the latter is sufficiently far from the plane boundary. The first is a solution very similar to Maxwell’s and the second pertains to an insulated sphere, in unbounded domains. We also introduce a simple correction to these approximations to account for the presence of the boundary, which improves their accuracy to smaller depths. These simple solutions provide insight into the physical mechanism on one hand and reduce the solution to two familiar cases for a wide range of the parameters values, on the other hand. The results of this paper may constitute the starting point for solving more general problems of heat conduction in nonhomogeneous media. A first generalization is for an arbitrary time dependence of the upper boundary temperature, which can be approached by Fourier or Laplace transforms. A second envisaged extension is for a media with multiple spherical inhomogeneities.

626

A. Rabinovich et al. / International Journal of Heat and Mass Transfer 55 (2012) 618–628

In turn, this may be instrumental in assessing the validity of the concept of effective thermal properties of the medium (see [11] Section 3.4.7). The approximations presented here may greatly simplify such a task, approached only by numerical methods in the past.

Appendix A. Derivation of the exact solution coefficients Substituting (14a,b) in (11a) yields

 1 X kn ðkex  r 0 Þ ðn þ 1=2ÞAn Pn ðlÞ þ ðn þ 1=2ÞAn ð1Þnþ1 kn ðkex Þ r¼1 n¼0 n¼0 1 X ~ Þ þ ekex ðldÞ ¼  P n ðl ðn þ 1=2ÞBn Pn ðlÞ ðA1Þ

1 X

n¼0

e ¼ p  l0 . Using the addition theorem [12] (see also [13]), where l

~Þ ¼ kn ðkr 0 ÞPn ðl

1 X 1 X

2p

i ð2s þ 1Þð2p þ 1Þ

Truncating the sums over n and s so that each series has N + 1 terms as opposed to infinity, allows writing Eq. (A7) in matrix form and solving a linear set of equations. With

0 B B B B B M¼B B B B B @

 where 1 X



j1 m1

n

s

p

0 0 0

j2 m2

j3 m3



ðn þ 1=2ÞAn Pn ðlÞ þ

2

ks ð2kdÞg p ðkrÞP p ðlÞ

ðA2Þ

is the Wigner 3j-symbol [14], we arrive at

1 X ¼ ðn þ 1=2ÞBn Pn ðlÞ

ðA3Þ

n¼0

 2 2p nþ1 n s p where Dnps ¼ i kð1Þ ð2s þ 1Þð2p þ 1Þ ks ð2kex dÞ. Applyn ðkex Þ 0 0 0 R1 ing 1 Pm ðlÞdl to all terms and taking advantage of Legendre orthogonality we obtain after some rearranging: 1 X 1 X 2n þ 1 Bm ¼ Am þ An g m ðkex ÞDnms þ Gm 2m þ1 n¼0 s¼0

ðA4Þ

þ lekex ðldÞ ¼

n¼0 p¼0 s¼0

rffiffiffiffi

~m ¼ 1 where k km ðkex Þ

1 jX

C h

ðn þ 1=2ÞBn g~n Pn ðlÞ

ðA5Þ

@km ðxÞ @x

x¼kex

and g~m ¼ g

1 m ðkin Þ

h

@g m ðxÞ @x

i x¼kin

pffiffiffiffi

jBm g~m 

pffiffiffiffi

jg~m ÞAm þ

ðA6Þ

pffiffiffiffi 1 X 1 g 0m ðkex Þ  jg~m  g m ðkex Þ X ð2n þ 1ÞAn Dnms ¼ Em 2m þ 1 n¼0 s¼0

pffiffiffiffi ffiffiffiffi b m ¼ R 1 ðpj g~m  lÞekex ðldÞ Pm ðlÞdl. The where Em ¼ jg~m Gm  G 1 b m have analytical expressions provided by [15]: integrals Gm and G 1

eðkex lÞ Pm ðlÞdl

1

¼ eðkex Þ

.. . 1 P e C N D0Ns

m h i X ðm þ lÞ! l 2 ðkex Þl1  ð1Þmþl  eð2kex Þ l!ðm  lÞ! l¼0

.. . 1 P e C N 3D1Ns

ð2N þ 1ÞDN1s

s¼0

.. . e N P ð2N þ 1ÞDNNs   C N þ C

s¼0

1

C C C C C C C C C C A

s¼0

ðA10Þ

EN

pffiffiffi ffiffiffiffi ~m  pj e m ¼ g0m ðkex Þ jg~m gm ðkex Þ. g~m and C where C m ¼ k 2mþ1

Appendix B. Derivation of the approximations through the exact solution

Neglecting the second term in (14a) results in the vanishing of the second term on the RHS of (A4) and (A6) to yields

Bm ¼ Am þ

Z

1

ekex ðldÞ Pm ðlÞdl

ðB1Þ

1

and

~m þ Am k

Z

1

pffiffiffiffi

lekex ðldÞ Pm ðlÞdl ¼ jBm g~m

ðB2Þ

1

We begin with the previous result of the large depth approximation given by (14a,b) and (16a,b). Then, we use the following relationships from [16]

p 2z

ez

n X ðn þ 1=2; kÞð2zÞk

ðB3Þ

0

where ðn; kÞ ¼ k!CCð1=2þnþkÞ and ð1=2þnkÞ

ðA7Þ

Z

e1 C

s¼0

1 0 1 E0 A0 BA C BE C B 1C B 1C C B C MB B .. C ¼ B .. C @ . A @ . A

g n ðzÞ ¼

b m ¼ R 1 lekex ðldÞ Pm ðlÞdl. Finally substituting (A4) in (A6) where G 1 and rearranging gives

~m  ðk

s¼0

s¼0 1 P

0

. After apply-

1 X 1 X 2n þ 1 bm An g 0m ðkex ÞDnms  G 2m þ1 n¼0 s¼0

1 e 1 P 3D11s   C1 þ C

the system becomes

kn ðzÞ ¼

ing Legendre orthogonality and rearranging:

~m ¼ Am k

s¼0

1

ðA9Þ

n¼0

i

e 0 P ð2N þ 1ÞDN0s C 1

 

B.2. The complex Maxwell approximation

1 X 1 X 1 X ~n Pn ðlÞ þ ðn þ 1=2ÞAn k ðn þ 1=2ÞAn g 0p ðkex ÞDnps Pp ðlÞ

n¼0

1

Solving this system of two equations for Am and Bm gives (16a,b).

R1 where Gm ¼ 1 ekex ðldÞ Pm ðlÞdl. Substituting (14a,b) in (11b) and utilizing (A2) we get 1 X

e 0 P 3D10s C

B.1. The large depth approximation

n¼0 p¼0 s¼0

þ ekex ðldÞ

s¼0 1 e 1 P D01s C

AN

1 X 1 X 1 X ðn þ 1=2ÞAn g p ðkex ÞDnps Pp ðlÞ

n¼0

1

s¼0

p¼0 s¼0



e 0 P D00s C0 þ C

ðA8Þ

zn 1  3  5    ð2n þ 1Þ " # z2 ðz2 Þ2 þ    ðB4Þ þ 2  1þ 2  1!ð2n þ 3Þ! 2  2!ð2n þ 3Þ!ð2n þ 5Þ!

to arrive at

Pn k kn ðkex  rÞ ekex ðr1Þ 0 ðn þ 1=2; kÞð2kex  rÞ ¼  P n k kn ðkex Þ r 0 ðn þ 1=2; kÞð2kex Þ

ðB5Þ

and

2

3 ðkin  rÞ2 ððkin  rÞ2 Þ2 þ   7 þ 61 þ 6 7 2  1!ð2n þ 3Þ 22  2!ð2n þ 3Þð2n þ 5Þ g n ðkin  rÞ 7 ¼ rn 6 6 7 2 2 2 g n ðkin Þ ðkin Þ ððkin Þ Þ 4 5 1þ þ  þ 2 2  1!ð2n þ 3Þ 2  2!ð2n þ 3Þð2n þ 5Þ ðB6Þ

627

A. Rabinovich et al. / International Journal of Heat and Mass Transfer 55 (2012) 618–628

Applying a perturbation expansion in (B5) for jkexj ? 0 and retaining only first order terms yields

kn ðkex  rÞ 1 ¼ nþ1 kn ðkex Þ r

ðB7Þ

Similarly, expanding (B6) for jkinj ? 0 and retaining only first order terms yields

g n ðkin  rÞ ¼ rn g n ðkin Þ

ðB8Þ

Next we expand (16a,b) in the same manner and use the following relationships given by [16]

@ nkn1 ðzÞ þ ðn þ 1Þknþ1 ðzÞ kn ðzÞ ¼  @z 2n þ 1

ðB9Þ

@ ng ðzÞ þ ðn þ 1Þg nþ1 ðzÞ g ðzÞ ¼ n1 @z n 2n þ 1

ðB10Þ

to arrive at

~n ¼  k

Pn1

k¼0 ðn

k

nþ1 2nþ1

Pnþ1

k

 1=2; kÞð2kex Þ þ k¼0 ðn þ 1:5; kÞð2kex Þ Pn k k¼0 ðn þ 1=2; kÞð2kex Þ

ðB11Þ and ðnþ2Þk2

g~n ¼

kin þ

k5in

ðB12Þ

þ 22 2!ð2nþ3Þð2nþ5Þ þ   

Expanding (B11) and (B12) for jkex/inj ? 0 and retaining only first order terms yields

~n ¼  ðn þ 1Þ k kex

ðB13Þ

and

g~n ¼

n ; kin

n > 0 and g~0 ¼

T 01;in  T 01;ex ¼ 1;

ðC3Þ

j

@T 01;in

T 02;in  T 02;ex ¼ kex lr;

@r

j



@T 01;ex

@T 02;in @r

@r 

¼ 0;

@T 02;ex @r

r¼1

¼ kex l;

ðC4Þ r¼1

ðC5Þ

Eqs. (C1) and (C2) are Laplace equations, thus having a known harmonic series solution given by [10]:

A0 A1 A2 þ 2 l þ 3 ð3l2  1Þ þ    ; r P 1 r r 2r 1 ¼ B0 þ B1 rl þ B2 r2 ð3l2  1Þ þ    ; r 6 1 2

T 01=2;ex ¼

ðC6Þ

T 01=2;in

ðC7Þ

kin 3

C.2. The insulated sphere approximation This time we apply a perturbation expansion to the external field only, given by (9a) and (10b), assuming jkexj ? 0, j ? 0 and retaining only first order terms to arrive at the external solution of (20a):

T 0IS;ex ¼ T 01;ex þ T 02;ex ¼

kex l ; rP1 2 r2

ðC8Þ

As for the internal field, we solve the Helmholtz equation (9a):

ðnþ4Þk4

in in n þ 21!ð2nþ3Þ þ 22 2!ð2nþ3Þð2nþ5Þ þ 

k3in 21!ð2nþ3Þ

r!1

Finally the solutions (C6) and (C7) with the conditions (C4) and (C5) gives rise to (18).

and

n 2nþ1

T 01;2 ¼ 0;

ðB14Þ

Finally substituting (B7) and (B8) in (14a,b) and (B13) and (B14) in (16a,b), expanding for jkexj ? 0 and retaining only first order terms, yields (18). B.3. The insulated sphere approximation We begin with (16a,b) by substituting (B13), applying j ? 0, expanding for jkexj ? 0 and retaining only first order terms to remain with A1 = kex/3, B0 = 2 and B1 = kex. Substituting these along with (B7) in (14a,b) yields (20a,b).

Appendix C. Derivation of the approximations by perturbation expansion of the general equations C.1. The complex Maxwell approximation The general equations to be solved far from the surface boundary are given by (9a,b), (10b) and (11a,b). If we apply a perturbation expansion in kex ; T 0 ¼ T 00 þ T 01 þ T 02 þ    where T 00 ¼ 0, to Eqs. (9a,b) and (10b), assuming jkex/inj ? 0 and retaining only first order terms we get:

r2 T 01 ¼ 0

ðC1Þ

r2 T 02 ¼ 0

ðC2Þ

r2 T 0in  k2in  T 0in ¼ 0;

z 6 d; r 6 1

ðC9Þ

with the following boundary conditions of given temperature from (C4) and (C5)

T 01;in ¼ 1 þ T 01;ex ¼ 1 ðr ¼ 1Þ T 02;in ¼ kex lr þ T 02;ex ¼

3 kex l ðr ¼ 1Þ 2

ðC10Þ ðC11Þ

and arrive at the solution of (20b). The work for this paper was carried out in partial fulfillment of the requirements for a doctor of philosophy degree by Avinoam Rabinovich at the Tel Aviv University. References [1] H. Carslaw, J. Jaeger, Conduction of Heat in Solids, second ed., Oxford University Press, London, 1959. [2] J. Maxwell, A treatise on electricity and magnetism, Oxford University Press (Clarendon), 1873. [3] E. Gordeliy, S. Crouch, S. Mogilevskaya, Transient heat conduction in a medium with two circular cavities: Semi-analytical solution, Int. J. Heat Mass Transfer 51 (2008) 3556–3570. [4] E. Gordeliy, S. Crouch, S. Mogilevskaya, Transient heat conduction in a medium with multiple spherical cavities, Int. J. Numer. Methods Eng. 77 (2009) 751– 775. [5] D. Oliver, Analytical solution for heat conduction near an encapsulated sphere with heat generation, J. Heat Transfer 130 (2008) 024502. [6] S. Traytak, On the time-dependent diffusive interaction between stationary sinks, Chem. Phys. Lett. 453 (2008) 212–216. [7] D. Shrivastava, R. Roemer, An analytical study of heat transfer in a finite tissue region with two blood vessels and general Dirichlet boundary conditions, Int. J. Heat Mass Transfer 48 (2005) 4090–4102. [8] I. Dincer, Estimation of dimensionless temperature distributions in spherical products during hydrocooling, Int. Commun. Heat Mass Transfer 22 (1995) 123–131. [9] J. Khedari, P. Benigni, J. Rogez, J.C. Mathieu, New apparatus for thermal diffusivity measurements of refractory solid materials by the periodic stationary method, Rev. Sci. Instrum. 66 (1995) 193–198. [10] G. Arfken, Mathematical Methods for Physicists, Academic Press, New York, London, 1970. [11] G. Dagan, Flow and Transport in Porous Formations, Springer-Verlag, Berlin, Heidelberg, 1989.

628

A. Rabinovich et al. / International Journal of Heat and Mass Transfer 55 (2012) 618–628

[12] H. Clercx, P. Schram, An alternative expression for the addition theorems of spherical wave solutions of the Helmholtz equation, J. Math. Phys. 34 (1993) 5292–5302. [13] T. Miloh, Hydrodynamics of deformable contiguous spherical shapes in an incompressible inviscid fluid, J. Eng. Math. 11 (1977) 349–372. [14] A. Messiah, Quantum Mechanics I, North-Holland, Amsterdam, 1981.

[15] Y.A. Cengel, M.N. Özisßik, Integrals involving Legendre polynomials that arise in the solution of radiation transfer, J. Quant. Spectrosc. Radiat. Transfer 31 (1984) 215–219. [16] M. Abramowitz, I. Stegun, Handbook of Mathematical Functions, Dover Publications, Inc., New York, 1965.