Two-dimensional Problems for Thermoelasticity, with Two Relaxation Times in Spherical Regions under Axisymmetric Distributions

PERGAMON

International Journal of Engineering Science 37 (1999) 299±314

Two-dimensional problems for thermoelasticity, with two relaxation times in spherical regions under axisymmetric distributions Hany H. Sherief *, Fouad A. Megahed Department of Mathematics, Faculty of Science, University of Qatar, Doha, Qatar Received 4 March 1998 (Communicated by E. S. S° UHUBIÇ )

Abstract Two-dimensional (2D) axisymmetric problems are considered within the context of the theory of thermoelasticity, with two relaxation times. The general solution is obtained in the Laplace transform domain by using a direct approach without the use of potential functions. The resulting formulation is utilized to solve a problem for a thick spherical shell. The surface of the shell is taken as traction free and subjected to given axisymmetric temperature distributions. The inversion of the Laplace transforms are carried out using the inversion formula of the transform together with Fourier expansion techniques. Numerical solutions are obtained for the temperature, displacement and stress distributions in the physical domain. Numerical results are represented graphically. # 1998 Elsevier Science Ltd. All rights reserved.

1. Introduction In 1967, the theory of generalized thermoelasticity with one relaxation time was introduced by Lord and Shulman [1]. The motivation behind the introduction of this theory was to deal with the apparent paradox of in®nite speeds of propagation predicted by the coupled theory of thermoelasticity, introduced by Biot [2] in 1956. The generalized equation of heat conduction is hyperbolic and, hence, automatically ensures ®nite speeds of wave propagation. This theory was extended [3] by Dhaliwal and Sherief to anisotropic media. Among the contributions to this theory are the proofs of uniqueness theorems by Ignaczak [4] and by Sherief [5]. The state space formulation for one-dimensional (1D) problems * Corresponding author. Tel.: 00-974-419134; Fax: 00-974-835061; E-mail: [email protected]. 0020-7225/98/$19.00 # 1998 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 2 2 5 ( 9 8 ) 0 0 0 7 0 - 6

300

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

was completed by Anwar and Sherief in Ref. [6] and by Sherief in Ref. [7]. The state space formulation for two-dimensional (2D) problems was carried out by Sherief and Anwar in Ref. [8]. The boundary element formulation was conducted by Anwar and Sherief in Ref. [9]. Sherief and Anwar have also solved a 2D problem of a thick plate with a moving heat source and a 2D problem for an in®nite cylinder, in Refs [10, 11], respectively. Sherief and Hamza have solved a 2D problem of a thick plate under axisymmetric temperature distribution, and discussed wave propagation for this theory in Ref. [12]. Green and Lindsay [13] developed the theory of thermoelasticity with two relaxation times, which is based on a generalized inequality of thermodynamics. This theory does not violate Fourier's law of heat conduction when the body under consideration has a center of symmetry. In this theory, both the equations of motion and of heat conduction are hyperbolic, but the equation of motion is modi®ed and di€ers from that of the coupled thermoelasticity theory. This theory was initiated by MuÈller [14]. It was further extended by Green and Laws [15]. The ®nal form used in the present work is that of Green and Lindsay [13]. This theory was also obtained independently by S° uhubi [16]. Longitudinal wave propagation for this theory was studied by Erbay and S° uhubi in Ref. [17]. Ignaczak proved a decomposition theorem in Ref. [18]. Sherief obtained the fundamental solution for this theory in Ref. [19], formulated the state space approach in Ref. [20], and solved a thermo-mechanical shock problem in Ref. [21]. The boundary integral equation formulation was carried out by Anwar and Sherief in Ref. [22]. Solutions of thermoelastic problems for spherical regions are not as numerous as those for Cartesian and cylindrical problems. Most of the treated problems are either 1D spherically symmetric ones or axisymmetric 2D problems under simplifying assumptions. Sternberg and Chakravorty [23] have solved a thermal shock uncoupled 1D problem. Hata has solved a coupled 1D thermal shock problem for a hollow sphere, caused by rapid uniform heating, in Ref. [24]. The general solution for spherically symmetric problems, with a heat source in generalized thermoelasticity valid for short times, was obtained by Sherief in Ref. [25]. Axially symmetric steady-state 2D problems in spherical regions were solved by McDowell and Sternberg in Ref. [26] and by Ignaczak in Ref. [27]. Ignaczak [28] and Piechocki [29] have solved dynamic problems in thermoelasticity by assuming that the time variable is harmonic, which tends to obscure the transient nature of the problems considered. Tanigawa and Kosako [30] have solved a transient coupled axially symmetric thermal stress problem for an in®nite medium with a spherical cavity, by neglecting inertia terms in their solution. Tanigawa and Takeuti [31] have obtained the three dimensional (3D) solution to coupled thermolelastic problems in spherical regions, again by neglecting inertia terms.

2. Formulation of the problem We consider a homogeneous isotropic thermoelastic solid under axisymmetric conditions. The axis of symmetry is the z-axis, and the origin of the system of coordinates is at the center of the sphere. We denote by (r, W, j) the spherical polar coordinates, and by t the time variable. Let T(r, W, t) be the axisymmetric temperature distribution throughout the solid. The

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

equations of motion in the absence of body forces have the form:   @ 2 ui @T;i ; r 2 ˆ …l 0 ‡ m 0 †uj;ij ‡ m 0 mi;jj ÿ g T;i ‡  @t @t

301

…1†

where u i = (u, v, 0) are the components of the displacement vector, r is the density, l 0 and m 0 are Lame's moduli,  is a constant that acts as a relaxation time, and g is a material constant given by g = (3l 0 + 2m 0 ) at, where at is the coecient of linear thermal expansion. The constitutive equations are given by:   @T 0 0 dij ; sij ˆ l ekk dij ‡ 2m eij ÿ l T ÿ T0 ‡  …2† @t where sij are the components of the stress tensor, eij are the components of the strain tensor and dij is Kronecker's delta symbol. The energy equation in the absence of heat source has the form: kT;ii ˆ rcE

@T @2 T @ui;i ‡ rtcE 2 ‡ gT0 ; @t @t @t

…3†

where k is the thermal conductivity, cE is the speci®c heat for processes with invariant strain tensor, t is a constant that acts as another relaxation time, and T0 is a reference temperature chosen such that v(T ÿ T0)/T0v W 1. Using the non-dimensional variables: r 0 ˆ c1 Zr;

u 0 ˆ c1 Zu;

 0 ˆ c1 Z;

t 0 ˆ c21 Zt;

sij …T ÿ T0 † ; yˆ ; 0 T0 m p where Z = rcE/k and c1 = …l 0 ‡ 2m 0 †=r is the speed of propagation of isothermal waves. Eqs. (1)±(3) takes the following form (dropping the primes for convenience):       2 1 @ @ @u @ @y 2@ u 2 @e ÿb ; b 2 ˆb ÿ sin W …rv† ÿ y‡ @t @r r2 sin W @W @r @W @r @t t 0 ˆ c21 Zt;

 0 ˆ c21 Z;

s 0ij ˆ

…4†

    @2 v b2 @e 1 @ @ @u @ @y b 2ˆ ‡ …rv† ÿ ÿ br y‡ ; @t @W @W @t r @W r @r @r

…5†

  @y dij ; sij ˆ …b2 ÿ 2†edij ‡ 2eij ÿ b y ‡  @t

…6†

2

@y @2 y @e ‡t 2 ‡g ; …7† @t @t @t where g = g/kZ, b 2 = (l 0 + 2m 0 )/m 0 , and b = gT0/m 0 . In the above equations e = u i,j is the r2 y ˆ

302

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

cubical dilatation given in spherical polar coordinates by: @u 2u 1 @…v sin W† ‡ ‡ ; …8† @r r r sin W @W and H 2 is Laplace's operator, given in spherical polar coordinates for the case of axisymmertic distributions by: eˆ

  @2 2 @ 1 @ @ ‡ sin W : r ˆ 2‡ r @r r2 sin W @W @W @r 2

Eqs. (4) and (5) can be combined to give:     @2 @y 2 2 r ÿ 2 e ˆ cr y ‡  ; @t @t

…9†

where c = b/b 2. We introduce the Laplace transform de®ned by the relation:  ˆ L‰f…t†Š ˆ f…s†

1 …

eÿst f…t†dt:

0

Assuming homogeneous initial conditions, then in the Laplace transform domain, Eqs. (7) and (9) can be written as:  …r2 ÿ s ÿ ts2 †y ˆ gse;

…10†

 …r2 ÿ s2 †e ˆ c…1 ‡ s†r2 y:

…11†

 Eliminating e between Eqs. (10) and (11), we obtain the following equation for y: fr4 ÿ ‰s2 ‡ s…1 ‡ ts† ‡ es…1 ‡ s†Šr2 ‡ s3 …1 ‡ ts†gy ˆ 0;

…12†

where e = cg is the coupling coecient. Eq. (12) can be factorized as: …r2 ÿ k21 †…r2 ÿ k22 †y ˆ 0;

…13†

where k 21 and k 22 are the roots of the characteristic equation: k4 ÿ ‰s2 ‡ s…1 ‡ ts† ‡ es…1 ‡ s†Šk2 ‡ s3 …1 ‡ ts† ˆ 0;

…14†

because of the linearity of Eq. (13), its solution is the sum of the solutions of the equations: …r2 ÿ k21 †y ˆ 0 and …r2 ÿ k22 †y ˆ 0: Solving the above two equations, we obtain Y in the form: y ˆ y 1 ‡ y 2 ; where

…15a†

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

303

1 2 X X y1 ˆ p1 Pn …m† …k2i ÿ s2 †Ani In‡12 …ki r†; r nˆ0 iˆ1

…15b†

1 2 X 1 X Pn …m† …k2i ÿ s2 †Bni Kn‡12 …ki r†: y 2 ˆ p r nˆ0 iˆ1

…15c†

In the above equations, In and Kn denote the modi®ed Bessel functions of the ®rst and second kinds of order n, respectively. Pn is the Legendre polynomial of degree n of argument m = cos W, and A ni, B ni are parameters depending on s only, that will be determined using the boundary conditions. The function y1 is bounded at the origin, while y2 is bounded at in®nity. Eliminating y between Eqs. (10) and (11), we obtain: …r2 ÿ k21 †…r2 ÿ k22 †e ˆ 0:

…16†

The solution of Eq. (16) that is compatible with Eq. (10), has the form: e ˆ e1 ‡ e2 ;

…17a†

where e1 ˆ

1 2 X c…1 ‡ s† X p Pn …m† k2i Ani In‡12 …ki r†; r nˆ0 iˆ1

1 2 X c…1 ‡ s† X e2 ˆ p Pn …m† k2i Bni Kn‡12 …ki r†: r nˆ0 iˆ1

…17b† …17c†

As before, e1 denotes the part of e that is bounded at the origin, while e2 denotes that part that is bounded at in®nity. In order to obtain the displacement component u, we shall use the Laplace transform of Eq. (4), namely: 2 @u 2u 2 @e @y ‡ 2 ÿ b2 s2 u ˆ e ‡ …1 ÿ b2 † ‡ b…1 ‡ s† : r @r r r @r @r Substituting for e and y from Eqs. (15) and (17) into Eq. (18), we obtain: r2 u ‡

r2 u ‡

1 2 @u 2u c…1 ‡ s† X ‡ 2 ÿ b2 s2 u ˆ Pn …m†‰f1 …r† ‡ f2 …r†Š: 3 r @r r r2 nˆ0

where f1 …r† ˆ

2 X iˆ1

f2 …r† ˆ

2 X iˆ1

Ani frki …k2i ÿ b2 s2 †In‡32 …ki r† ‡ ‰…n ‡ 2†k2i ÿ nb2 s2 ŠIn‡12 …ki r†g; Bni fÿrki …k2i ÿ b2 s2 †Kn‡32 …ki r† ‡ ‰…n ‡ 2†k2i ÿ nb2 s2 ŠKn‡12 …ki r†g:

…18†

…19†

304

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

In obtaining the above expressions, we have used the following properties of the modifed Bessel functions [32]: d n ‡ 1=2 In‡1=2 …kr† ˆkIn‡3=2 …kr† ‡ In‡1=2 …kr†; dr r d n ‡ 1=2 Kn‡1=2 …kr† ˆ ÿkKn‡3=2 …kr† ‡ Kn‡1=2 …kr†: dr r

…20b†

The solution of Eq. (19) can be written as: u ˆ u1 ‡ u2 ;

…21a†

where u1 ˆ

u2 ˆ

( 1 c…1 ‡ s† X 3

r2

Pn …m†

2 X

nˆ0

iˆ1

( 1 c…1 ‡ s† X

2 X

3

r2

nˆ0

Pn …m†

iˆ1

Ani ‰rki In‡32 …ki r† ‡ nIn‡12 …ki r†Š ‡

1 X nˆ1

Bni ‰ÿrki Kn‡32 …ki r† ‡ nKn‡12 …ki r†Š ‡

) Cn Pn …m†In‡12 …bsr† ;

1 X nˆ1

…21b† )

Dn Pn …m†Kn‡12 …bsr† ; …21c†

where C n and D n are parameters depending on s only. In obtaining the above solutions we have made use of the following properties of the modi®ed Bessel functions [32]: d n ‡ 3=2 In‡3=2 …kr† ˆKIn‡1=2 …kr† ÿ In‡3=2 …kr†; dr r d n ‡ 3=2 Kn‡3=2 …kr† ˆ ÿ kKn‡1=2 …kr† ÿ Kn‡3=2 …kr†: dr r

…22a†

…22b†

In order to obtain the displacement component v, we shall use the Laplace transform of Eq. (8), for the cubical dilatation, namely: e ˆ

1 @ 2 1 @  ‡ …r u† …v sin W†; 2 r @r r sin W @W

…23†

which implies that:   @ @u 2 …v sin W† ˆ r ‡ u ÿ e : @m @r r

…24†

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

305

Substituting on the right-hand side of Eq. (24) for u and e from Eqs. (21) and (17), then integrating the resulting euation with respect to m, we obtain: v ˆ v1 ‡ v2 ;

…25a†

where (  ) 1 2 X c…1 ‡ s† X bsr I 3 …bsr† ; v1 ˆ 3 ‰mPn …m† ÿ Pnÿ1 …m†Š nAni In‡12 …ki r† ‡ Cn In‡12 …bsr† ‡ n ‡ 1 n‡2 r2 sin W nˆ1 iˆ1 …25b† (  ) 1 2 X c…1 ‡ s† X bsr K 3 …bsr† : ‰mPn …m† ÿ Pnÿ1 …m†Š nBni Kn‡12 …ki r† ‡ Dn Kn‡12 …bsr† ÿ v2 ˆ 3 n ‡ 1 n‡2 r2 sin W nˆ1 iˆ1 …25c† We have used the following integral relation of the Legendre polynomials [32] in obtaining Eq. (25)a)±(25c) … mPn …m† ÿ Pnÿ1 …m† : …26† Pn …m†dm ˆ n‡1 It can be easily shown that v is bounded at W = 0, p, 2p. In fact: lim v ˆ lim v ˆ lim v ˆ 0:

W40

W4p

W42p

From Eq. (2), it follows that the stress components in the Laplace transform domain have the following form: @u  ‡ …b2 ÿ 2†e ÿ b…1 ‡ s†y; @r 2 @v 2u  ‡ ‡ …b2 ÿ 2†e ÿ b…1 ‡ s†y; s WW ˆ r @W r 2 cot W 2u  s ff ˆ v ‡ ‡ …b2 ÿ 2†e ÿ b…1 ‡ s†y; r r 1 @u v @v ÿ ‡ ; s rW ˆ r @W r @r s rf ˆs Wf ˆ 0: s rr ˆ2

…27a† …27b† …27c† …27d† …27e†

Substituting from Eqs. (15), (17) and (21), into Eq. (27)a), we obtain: s rr ˆ s rr1 ‡ s rr2 ; where

…28a†

306

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

s rr1 ˆ

 1 2c…1 ‡ s† X 5

r2 ‡

1 X nˆ1

Pn …m†

nˆ0

2 X

 Ani

iˆ1

  1 n…n ÿ 1† ‡ b2 s2 r2 In‡12 …ki r† ÿ 2rki In‡32 …ki r† 2

 Cn Pn …m†‰…n ÿ 1†In‡12 …bsr† ‡ bsrIn‡32 …bsr†Š ;

…28b†

and s rr2 ˆ

 1 2c…1 ÿ s† X 5

r2 ‡

1 X nˆ1

nˆ0

Pn …m†

2 X

 Bni

iˆ1

  1 2 2 2 n…n ÿ 1† ‡ b s r Kn‡12 …ki r† 2

 Dn Pn …m†‰…n ÿ 1†Kn‡12 …bsr† ÿ bsrKn‡32 …bsr†Š :

…28c†

Substituting from Eqs. (21) and (25) into Eq. (27)d), we obtain: s rW ˆ s rW1 ‡ s rW2 ;

…29a†

where s rW1 ˆ

s rW2

X 1 2 2c…1 ‡ s† X ‰mP …m† ÿ P …m†Š nAni ‰…n ÿ 1†In‡12 …ki r† ‡ ki rIn‡32 …ki r†Š n nÿ1 5 r2 sin W nˆ1 iˆ1   Cn 1 2 2 2 2 …n ÿ 1 ‡ b s r †In‡12 …bsr† ÿ bsrIn‡32 …bsr† ; ‡ 2 …n ‡ 1†

X 1 2 2c…1 ‡ s† X ˆ 5 ‰mPn …m† ÿ Pnÿ1 …m† nBni ‰…n ÿ 1†Kn‡12 …ki r† ÿ ki rKn‡32 …ki r†Š r2 sin W nˆ1 iˆ1   Dn 1 ‡ …n2 ÿ 1 ‡ b2 s2 r2 †Kn‡12 …bsr†Kn‡32 …bsr† ; 2 …n ‡ 1†

…29b†

…29c†

where we have used the relation [32]: dPn …cos W† n ˆ ‰cos WPn …cos W† ÿ Pnÿ1 …cos W†Š: dW sin W

…30†

The remaining stress components can be obtained in a similar manner. Eqs. (15), (21), (25), (28) and (29) give the general solution of the problem in the Laplace transform domain in terms of the parameters A ni, B ni, C n and D n. These parameters can be obtained from the boundary condition of the particular problem under consderation.

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

307

3. Application In order to illustrate the above results, we shall solve a problem for a thick spherical shell whose surfaces are spheres with centers at the origin and of radii a and b (a < b). The surfaces of the shell are assumed to be traction free. The inner and outer surfaces of the shell are in contact with media of known axisymmetric temperatures F1(m, t) and F2(m, t), with heat transfer coecients of L1 and L2, respectively. The conditions that the inner and outer surfaces are traction free, give, in the Laplace transform domain, the four equations: s rr …a; W; s† ˆ0;

…31†

s rr …b; W; s† ˆ0;

…32†

s rW …a; W; s† ˆ0;

…33†

s rW …b; W; s† ˆ0:

…34†

The thermal boundary conditions, give, on the inner and outer surfaces, respectively: qr …a; W; t† ˆL1 ‰F1 …m; t† ÿ y…a; W; t†Š;

…35†

qr …b; W; t† ˆL2 ‰y…a; W; t† ÿ F2 …m; t†Š;

…36†

where qr is the heat ¯ux in the radial direction. Using Fourier's law of heat conduction, which is valid for the theory of thermoelasticity with two relaxation times [13], namely: qn ˆ ÿk

@T ; @n

where qn is the component of the heat ¯ux vector in direction of the outward normal n. Using the non-dimensional variables used earlier, together with the non-dimensional heat ¯ux vector q 0 n = qn/kT0c1Z and dropping the prime, we obtain: qn ˆ ÿ

@y : @n

Eqs. (35) and (36) in the Laplace transform domain, thus, take the following form:  W; s† @y…a;  W; s† ˆ ÿL1 F1 …m; s†; ÿ L1 y…a; @r

…37†

 W; s† @y…b;  W; s† ˆ L2 F2 …m; s†: ‡ L2 y…b; @r

…38†

308

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

Expanding the functions F1(m, s) and F2(m, s) in a series of Legendre polynomials, we obtain: Fi …m; s† ˆ

1 X

fni …s†Pn …m†;

i ˆ 1; 2;

…39a†

nˆ0

where 2n ‡ 1 fni …s† ˆ 2

…1

Fi …m; s†Pn …m†dm;

i ˆ 1; 2; n ˆ 0; 1; 2; . . .

…39b†

ÿ1

Di€erentiating both sides of the equations and using Eq. (20a) and (20b), we obtain: 1 2 nn o X @y 1 1 X Pn …m† Ani …k2i ÿ s2 † In‡12 …ki r† ‡ ki In‡32 …ki r† ; ˆ p r @r r nˆ0 iˆ1

…40a†

1 2 nn o X @y 2 1 X Pn …m† Bni …k2i ÿ s2 † Kn‡12 …ki r† ÿ ki Kn‡32 …ki r† : ˆ p r @r r nˆ0 iˆ1

…40b†

Using Eqs. (28), (40) and (15), Eqs. (31), (32), (37) and (38) give upon equating coecients of P0(m) on both sides: 2 X p …k2i ÿ s2 †fA0i ‰L1 I12 …ki a† ÿ ki I32 …ki a†Š‡B0i ‰L1 K12 …ki a† ‡ ki K32 …ki a†Šg ˆ L1 af01 …s†;

…41†

2 X p …k2i ÿ s2 †fA0i ‰L2 I12 …ki b† ‡ ki I32 …ki b†Š‡B0i ‰L2 K12 …ki b† ÿ ki K32 …ki b†Šg ˆ L2 bf02 …s†;

…42†

iˆ1

iˆ1

2 X fA0i ‰b2 s2 a2 I12 …ki a† ÿ 4aki I32 …ki a†Š‡B0i ‰b2 s2 a2 K12 …ki a† ‡ 4aki K32 …ki a†Šg ˆ 0;

…43†

2 X fA0i ‰b2 s2 b2 I12 …ki b† ÿ 4bki I32 …ki b†Š‡B0i ‰b2 s2 b2 K12 …ki b† ‡ 4bki K32 …ki b†Šg ˆ 0:

…44†

iˆ1

iˆ1

Eqs. (41)±(44) are a set of four linear equations in the four unknown parameters A 01, A 02, B 01 and B 02, which can be easily solved to give the four unknowns. We will omit the details here.

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

309

For n > 0, Eqs. (15), (28), (29) and (40), together with Eqs. (31)±(34), (37) and (38), give upon equating coecients of P n(m) on both sides: 2 X iˆ1

   1 2 2 2 Ani n…n ÿ 1† ‡ b s a In‡12 …ki a† ÿ 2aki In‡32 …ki a† ‡ Cn ‰…n ÿ 1†In‡12 …bsa† ‡ bsaIn‡32 …bsa†Š 2

‡

2 X

 Bni

iˆ1

  1 n…n ÿ 1† ‡ b2 s2 a2 Kn‡12 …ki a† ‡ 2aki Kn‡32 …ki a† 2

‡ Dn ‰…n ÿ 1†Kn‡12 …bsa† ÿ bsaKn‡32 …bsa†Š ˆ 0; 2 X iˆ1

…45†

   1 2 2 2 Ani n…n ÿ 1† ‡ b s b In‡12 …ki b† ÿ 2bki In‡32 …ki b† ‡ Cn ‰…n ÿ 1†In‡12 …bsb† ‡ bsbIn‡32 …bsb†Š 2

‡

2 X

 Bni

iˆ1

  1 2 2 2 n…n ÿ 1† ‡ b s b Kn‡12 …ki b† ‡ 2bki Kn‡32 …ki b† 2

‡ Dn ‰…n ÿ 1†Kn‡12 …bsb† ÿ bsbKn‡32 …bsb†Š ˆ 0; 2 X iˆ1

h i nAni …n ÿ 1†In‡12 …ki a† ‡ ki aIn‡32 ‡

‡

2 X iˆ1

…46†

  Cn 1 …n2 ÿ 1 ‡ b2 s2 a2 †In‡12 …bsa† ÿ bsaIn‡32 …bsa† 2 …n ‡ 1†

nBni ‰…n ÿ 1†Kn‡12 …ki a† ÿ ki aKn‡32 …ki a†Š

  Dn 1 2 2 2 2 ‡ n ÿ 1 ‡ b s a †Kn‡12 …bsa† ‡ bsaKn‡32 …bsa† ˆ 0; 2 …n ‡ 1† 2 X iˆ1

Cn nAni ‰…n ÿ 1†In‡12 …ki b†In‡32 …ki b†Š ‡ …n ‡ 1†

‡

2 X iˆ1



…47†

  1 2 2 2 n ÿ 1 ‡ b s b In‡12 …bsb† ÿ bsbIn‡32 …bsb† 2 2

nBni ‰…n ÿ 1†Kn‡12 …ki b† ÿ ki bKn‡32 …ki b†Š

Dn ‡ …n ‡ 1†



  1 2 2 2 n ÿ 1 ‡ b s b Kn‡12 …bsb† ‡ bsbKn‡32 …bsb† ˆ 0; 2 2

…48†

310

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

 h 2 i hn i   X n 2 2 ÿ L1 In‡12 …ki a† ‡ ki In‡32 …ki a† ‡ Bni ÿ L1 Kn‡12 …ki a† ÿ ki Kn‡32 …ki a† …ki ÿ s † Ani a a iˆ1 p ˆ ÿL1 afn1 …s†;

…49†

2 n hn i hn io   X ‡ L2 In‡12 …ki b† ‡ ki In‡32 …ki b† ‡ Bni ‡ L2 Kn‡12 …ki b† ÿ ki Kn‡32 …ki b† …k2i ÿ s2 † Ani b b iˆ1

p ˆ L2 bfn2 …s†:

…50†

Eqs. (45)±(50) are a set of six linear equations in the six unknown parameters A n1, A n2, B n1, B n2, C n and D n, which can be solved to give the six unknowns. We will omit the details here. Actually, both linear systems are solved numerically during the process of the numerical inversion of the Laplace transform. This completes the solution of the problem in the Laplace transform domain.

4. Inversion of the Laplace transforms We shall now outline the numerical inversion method used to ®nd the solution in the physical domain. Let f(r, W, s) be the Laplace transform of a function f(r, W, t). The complex inversion formula for Laplace transforms can be written as: 1 f…r; W; t† ˆ 2pi

d‡i1 …

 W; s†ds; est f…r;

dÿi1

where d is an arbitrary real number greater than all the real parts of the singularities of f(r, W, s). Taking s = d + iy, the above integral takes the form: edt f…r; W; t† ˆ 2p

1 …

 W; d ‡ iy†y: eity f…r;

ÿ1

Expanding the function h(r, W, t) = exp(ÿdt)f(r, W, t) in a Fourier series in the interval [0, 2 T], we obtain the approximate formula [33]: f…r; W; t† ˆ f1 …r; W; t† ‡ ED ; where 1 X 1 f1 …r; W; t† ˆ c0 …r; W; t† ‡ ck …r; W; t†; 2 kˆ1

for

0RtR2T;

…51†

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

311

and ck …r; W; t† ˆ

edt  W; d ‡ ikp=T†Š; Re‰eikpt=T f…r; T

k ˆ 0; 1; 2

…52†

ED, the discretization error, can be made arbitrarily small by choosing the constant d large enough [33]. Since the in®nite series in Eq. (51) can only be summed up to a ®nite number N of terms, the approximate value of f(r, W, t) becomes: N X 1 fN …r; W; t† ˆ c0 …r; W; t† ‡ ck …r; W; t†; 2 kˆ1

for

0RtR2T:

…53†

Using the above formula to evaluate f(r, W, t), we introduce a truncation error ET that must be added to the discretization error to produce the total approximation error. Two methods are used to reduce the total error. First, the ``Korrecktur'' method [33] is used to reduce the discretization error. Next, the e-algorithm is used to reduce the truncation error and, hence, to accelerate convergence. The Korrecktur method uses the following formula to evaluate the function f(r, W, t): f…r; W; t† ˆ f1 …r; W; t† ÿ e2dT f1 …r; W; 2T ‡ t† ‡ E 0D ; where the discretization error vE 0 DvWvEDv [33]. Thus, the approximate value of f(r, W, t) becomes: fNK …r; W; t† ˆ fN …r; W; t† ÿ eÿ2dT fN 0 …r; W; 2T ‡ t†;

…54†

N 0 is an integer such that N 0 < N. We shall now describe the e-algorithm that is used to accelerate the convergence of the series in Eq. (53). Let N be an odd natural number, and let: sm …r; W; t† ˆ

m X

ck …r; W; t†

kˆ1

be the sequence of partial sums of Eq. (53). We de®ne the E-sequence by: e0;m ˆ 0; e1;m ˆ sm

Table 1 The constants of the problem k = 386 m 0 = 3.86  1010 b2=4 E = 0.0168 L1 = 1

at = 1.78  10 ÿ 5 l 0 = 7.76  1010 T0 = 293 t = 0.02 L2 = 2

cE = 383.1 r = 8954 b = 0.042  = 0.02 a= 1

Z = 8886.73 c1 = 4.158  103 g = 1.61 c = 0.01 b= 2

312

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

Fig. 1. Temperature distribution.

and ep‡1;m ˆ epÿ1;m‡1 ‡ 1=…ep;m‡1 ÿ ep;m †;

p ˆ 1; 2; 3; . . .

It can be shown that [33] the sequence e1;1 ; e3;1 ; e5;1 ; . . . ; eN;1 converges to f(r, W, t) + ED ÿc0/2 faster than the sequence of partial sums sm ; m ˆ 1; 2; 3; . . . The actual procedure used to invert the Laplace transforms consists of using Eq. (54) together with the e-algorithm. The values of d and T are chosen according to the criteria outlined in Ref. [33].

Fig. 2. Displacement distribution.

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

313

Fig. 3. Radial stress distribution.

5. Numerical Results The copper material was chosen for purposes of numerical evaluations. The constants of the problem are shown in Table 1. The normalized temperatures of the inner and outer surroundings were taken, respectively, as: Fi …m; t† ˆ m2 ; F2 …m; t† ˆ m3 : We thus obtain from Eq. (39)a) and (39b): f01 …s† ˆ1=3s; f21 …s† ˆ 2=s; f11 …s† ˆ f31 …s† ˆ f41 …s† ˆ ˆ 0: f12 …s† ˆ3=5s; f32 …s† ˆ 2=5s; f02 …s† ˆ f22 …s† ˆ f12 …s† ˆ ˆ 0: The in®nite series in all the above solutions thus terminate after n = 3. The computations were performed for three values of non-dimensional time, namely t = 0.05, t = 0.1 and t = 0.2. The numerical technique outlined above was used to obtain the temperature, the displacement, and the radial stress distributions as functions of r for W = 0. In all the ®gures, solid lines represent the fuction when t = 0.05, lines with an asterisk marker represent the function when t = 0.1, while lines with a dot marker represent the function when t = 0.2. The normalized temperature increment, y, is represented by the graph in Fig. 1. The displacement u is shown in Fig. 2. The radial stress component s n is shown in Fig. 3. The phenomenon of ®nite speeds of propagation is manifested in all these ®gures. For the smallest values of time considered we see that the heat e€ects of the surrounding media is localized in a region adjacent to the walls. This region expands with the passage of time to ®ll the whole cylinder for large values of time. This region corresponds to the propagation of wave fronts from the surfaces of the shell. This is not the case when using the coupled equation of heat conduction, where the thermal e€ects extend to ®ll the whole shell immediately.

314

H.H. Sherief, F.A. Megahed / International Journal of Engineering Science 37 (1999) 299±314

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

H. Lord, Y. Shulman, J. Mech. Phys. Solid 15 (1967) 299±309. M. Biot, J. Appl. Phys. 27 (1956) 240±253. R.S. Dhaliwal, H. Sherief, Q. Appl. Math. 38 1980 1±8. J. Ignaczak, J. Thermal Stresses 5 (1982) 257±263. H. Sherief, Q. Appl. Math. 45 (1987) 773±778. M. Anwar, H. Sherief, J. Thermal Stresses 11 (1988) 353±365. H. Sherief, J. Thermal Stresses 16 (1993) 163±180. H. Sherief, M. Anwar, J. Thermal Stresses 17 (1994) 567±590. M. Anwar, H. Sherief, Appl. Math. Model 12 (1988) 161±166. H. Sherief, M. Anwar, J. Thermal Stresses 15 (1992) 489±505. H. Sherief, M. Anwar, J. Thermal Stresses 17 (1994) 213±227. H. Sherief, F. Hamza, J. Thermal Stresses 17 (1994) 435±452. A.E. Green, K.A. Lindsay, J. Elast 2 (1972) 1±7. I. MuÈller, Arch. Rat. Mech. Anal. 41 (1971) 319±332. A.E. Green, N. Laws, Arch. Rat. mech. Anal. 45 (1972) 47±53. E.S. S° uhubi, Thermoelastic solids. in: A.C. Eringen (Eds.), Continuum Physics, chap. 2, vol. 2, Academic, New York, 1975. S. Erbay, E.S. S° uhubi, J. Thermal Stresses 9 (1986) 279±295. J. Ignaczak, J. Thermal Stresses 1 (1978) 41±52. H. Sherief, Int. J. Engng Sci. 30 (1992) 861±870. H. Sherief, Int. J. Engng Sci. 31 (1993) 1177±1189. H. Sherief, Int. J. Engng Sci. 32 (1994) 313±325. E. Sternberg, J. Chakravorty, Q. Appl. Math. 16 (1959) 205±218. T. Hata, J. Appl. Mech. 58 (1991) 63±69. H. Sherief, J. Thermal Stresses 9 (1986) 151±164. E. McDowell, E. Sternberg, J. Appl. Mech. 24 (1957) 376±380. J. Ignaczak, Arch. Mech. Stos. 12 (1960) 415±435. J. Ignaczak, Arch. Mech. Stos. 11 (1959) 399±408. W. Piechocki, Arch. Mech. Stos. 12 (1960) 553±561. Y. Tanigawa, M. Kosako, Int. J. Engng Sci. 24 (1986) 309±321. Y. Tanigawa, Y. Takeuti, ZAMM 63 (1983) 317±324. I.S. Sokolniko€, Mathematical Theory of Elasticity. 2nd ed., McGraw-Hill, New York, 1956. M. Spiegel, Mathematical Handbook. McGraw-Hill, New York, 1996. G. Honig, U. Hirdes, J. Comp. Appl. Math. 10 (1984) 113±132.