- Email: [email protected]
A two-dimensional thermoelasticity problem for thermomechanical shock with two relaxation times
Applied Mathematics and Computation 170 (2005) 172–184 www.elsevier.com/locate/amc
A two-dimensional thermoelasticity problem for thermomechanical shock with two relaxation times Nasser M. El-Maghraby *, Hamdy M. Yossef Department of Basic and Applied Science, Arab Academy for Science and Technology, P.O. Box 1029, Alexandria, Egypt
Abstract In this work a two-dimensional problem of thermoelasticity with two relaxation times is introduced. Laplace and Fourier transform techniques are used. The resulting formulation is applied to a thermomechanical shock half-space problem. The solution in the transformed domain obtained by using a direct approach. Numerical inversion of both transforms is carried out to obtain the temperature, stress and displacement distributions in the physical domain. Numerical results are represented graphically. Ó 2004 Elsevier Inc. All rights reserved. Keywords: Elasticity; Thermoelasticity; Two relaxation times; Mechanical shock
1. Introduction Danilovskaya [1] was the first to solve an actual problem in the theory of elasticity with uniform heat. The problem was for a half-space subjected to a *
Corresponding author. E-mail address: [email protected] (N.M. El-Maghraby).
0096-3003/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.10.077
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
173
thermal shock in the context of what became known as the theory of uncoupled thermoelasticity. In this theory the temperature is governed by a parabolic partial differential equation which does not contain any elastic terms. It was not much latter that many attempts were made to remedy the shortcomings of this theory. Biot [2] formulated the theory of coupled thermoelasticity to eliminate the paradox inherent in the classical uncoupled theory that elastic changes have no effect on the temperature. The heat equations for both theories, however, are of the diffusion type predicting infinite speeds of propagation for heat waves contrary to physical observations. Lord and Shulman [3] introduced the theory of generalized thermoelasticity with one relaxation time for the special case of an isotropic body. This theory was extended by Dhaliwal and Sherief [4] to include the anisotropic case. In this theory a modified law of heat conduction including both the heat flux and its time derivative replaces the conventional FourierÕs law. The heat equation associated with this theory is hyperbolic and hence eliminates the paradox of infinite speeds of propagation inherent in both the uncoupled and coupled theories of thermoelasticity. Uniqueness of solution for this theory was proved under different conditions by Ignaczak [5,6], Sherief and Dhaliwal [7], Dhaliwal and Sherief [4] and Sherief [8]. Green and Lindsay [9] obtained the theory of generalized thermoelasticity with two relaxation times. The uniqueness of solution for this theory was proved by Green [10]. The fundamental solution was obtained by Sherief [11]. Sherief and Helmy have solved a two-dimensional problem in [12]. Sherief and El-Maghraby have solved a crack problem in one relaxation time [13], ElMaghraby [14] has solved a two-dimensional problem for a half space subjected to heat sources. The authors have solved a one-dimensional problem with thermomechanical shock [15]. In this work a two-dimensional problem of thermoelasticity with two relaxation times is introduced. Laplace and Fourier transform techniques are used. The resulting formulation is applied to a thermomechanical shock half-space problem. The solution in the transformed domain obtained by using a direct approach. Numerical inversion of both transforms is carried out to obtain the temperature, stress and displacement distributions in the physical domain. Numerical results are represented graphically.
2. Formulation of the problem We shall consider a homogeneous isotropic thermoelastic solid occupying the half-space y P 0. The y-axis is taken perpendicular to the bounding plane pointing inwards. We shall also assume that the initial state of the medium is quiescent. The surface of this medium is subjected to a thermomechanical shock. The displacement vector thus, has components.
174
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
u ¼ ðu; v; 0Þ and the cubical dilatation e is given by e ¼ div u ¼
ou ov þ : ox oy
ð1Þ
The governing equations for thermoelasticity with two relaxation times consist of: (1) The equation of motion in vector form o2 u oT q 2 ¼ ðk þ lÞ grad e þ lr2 u c grad T þ m ; ð2Þ ot ot which has two Cartesian components o2 u oe o2 u o2 u oT o2 T þm ; q 2 ¼ ðk þ lÞ þ l 2 þ l 2 c ot ox ox oy ox ox ot
ð3Þ
o2 v oe o2 v o2 v oT o2 T þ l þ m ¼ ðk þ lÞ þ l c ; ot2 oy ox2 oy 2 oy oy ot
ð4Þ
q
where k, l are Lame´Õs constants, q is the density, t is the time variable, T is the absolute temperature of the medium, m is a constant with dimensions of time, called a relaxation time, and c is a material constant given by c = (3k + 2l)at, at being the coefficients of linear thermal expansion. (2) The equation of heat conduction oT o2 T oe þ s 2 þ cT 0 ; ð5Þ kr2 T ¼ qC E ot ot ot where k is the thermal conductivity of the medium, T0 is a reference temperature, CE is the specific heat at constant strain and s is another relaxation time. Applying the div operator to both sides of Eq. (2), we obtain o2 e oT q 2 ¼ ðk þ 2lÞr2 e cr2 T þ m ; ð6Þ ot ot where $2 is LaplaceÕs operator. The above equations are supplemented with the constitutive equations ov oT rxx ¼ ðk þ 2lÞe 2l c T T 0 þ m ; ð7aÞ oy ot ou oT ryy ¼ ðk þ 2lÞe 2l c T T 0 þ m ; ð7bÞ ox ot
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
oT rzz ¼ ke c T T 0 þ m ; ot ou ov þ rxx ¼ l ; oy ox
175
ð7cÞ
ð7dÞ
rxz ¼ ryx ¼ 0:
ð7eÞ
We shall use the following non-dimensional variables x0 ¼ c1 gx;
y 0 ¼ c1 gy;
t0 ¼ c21 gt;
u0 ¼ c1 gu;
v0 ¼ c1 gv;
cðT T 0 Þ rij ; rij ¼ ; h¼ k þ 2l l pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where g = (qcE/k) and c1 ¼ ðk þ 2l=qÞ is the speed of propagation of isothermal elastic waves. Eq. (1) retains its form, while Eqs. (3)–(6) take the following form (dropping the primes for convenience). 2 oe o2 h 2o u 2 2 oh 2 þm b 2 ¼ ðb 1Þ þ r u b ; ð8Þ ot ox ox ox ot s0 ¼ c21 gs;
b2
m0 ¼ c21 gm;
o2 v oe o2 h 2 2 oh 2 þ r þ m ¼ ðb 1Þ v b ; ot2 oy oy ox ot
oh o2 h oe þs 2 þe ; ot ot ot o2 e oh b 2 2 ¼ r 2 e b2 r 2 h þ m ; ot ot r2 h ¼
ð9Þ
ð10Þ ð11Þ
where b2 ¼
k þ 2l l
and
e¼
T 0 c2 : qcE ðk þ 2lÞ
The constitutive equation (7) in non-dimensional form can be written as ov oh rxx ¼ b2 e 2 b2 h þ m ; ð12aÞ oy ot ou oh ryy ¼ b2 e 2 b2 h þ m ; ð12bÞ ox ot oh rzz ¼ ðb2 2Þe b2 h þ m ; ð12cÞ ot
176
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
rxy ¼
ou ov þ : oy ox
ð12dÞ
The above equations are solved subject to the initial conditions h ¼ u ¼ v ¼ 0;
oh ou ov ¼ ¼ ¼ 0 at t ¼ 0 ot ot ot
ð13Þ
and the boundary conditions rxy ¼ 0;
ryy ¼ f ðx; tÞ
h ¼ Gðx; tÞ
at y ¼ 0;
at y ¼ 0;
ð14Þ ð15Þ
where f(x, t) and G(x, t) are known functions. Using FourierÕs law of heat conduction in non-dimensional form (which is valid for Green theory), namely qn ¼
oh ; on
ð16Þ
where qn denotes the normal component of the heat flux vector.
3. Solution in the transformed domain Taking the Laplace transform with parameter s (denoted by a bar) of both sides of Eqs. (8)–(12), we obtain the following set of equations ðr2 s ss2 Þ h ¼ ese;
ð17Þ
ðr2 sÞe ¼ ð1 þ msÞr2 h;
ð18Þ
b 2 s2 u ¼ ðb2 1Þ
oe oh þ r2 u b2 ð1 þ msÞ ; ox ox
ð19Þ
b2 s2v ¼ ðb2 1Þ
oe o h þ r2v b2 ð1 þ msÞ ; oy oy
ð20Þ
xx ¼ b2e 2 r
ov b2 ð1 þ msÞ h; oy
ð21aÞ
yy ¼ b2e 2 r
o u b2 ð1 þ msÞ h; ox
ð21bÞ
xy ¼ r
o u ov þ oy ox
ð21cÞ
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
177
and the Laplace transform of the cubical dilation (1) becomes o u ov e ¼ þ : ð22Þ oy ox We have assumed that the Laplace transforms of the temperature, stress components, displacements, and the applied thermomechanical shock exist. We now introduce the exponential Fourier transform (denoted by an asterisk) with respect to the space variable x, defined by Z 1 1 f ðq; y; sÞ ¼ I½f ðx; y; sÞ ¼ pffiffiffiffiffiffi eiqx f ðx; y; sÞ dx 2p 1 with corresponding inversion formula Z 1 1 eiqx f ðq; y; sÞ dq f ðx; y; sÞ ¼ I1 ½f ðq; y; sÞ ¼ pffiffiffiffiffiffi 2p 1
where i ¼
pffiffiffiffiffiffiffi 1:
We assume that all the relevant functions (temperature, stress, . . ., etc.) are sufficient smooth on the real line such that the Fourier transforms of these functions exists. Taking the Fourier transform of both sides of Eqs. (17)–(22), we obtain the following set of equations ðD2 q2 s ss2 Þ h ¼ ese ; ð23Þ ðD2 q2 s2 Þe ¼ ð1 þ msÞðD2 q2 Þ h;
ð24Þ
u ¼ iq½b2 ð1 þ msÞ h ðb2 1Þe ; ðD2 a2 Þ
ð25Þ
h ðb2 1ÞDe ; ðD2 a2 Þv ¼ b2 ð1 þ msÞD
ð26Þ
xx ¼ b2e 2Dv b2 ð1 þ msÞ h; r
ð27aÞ
yy ¼ b2e 2iq r u b2 ð1 þ msÞ h;
ð27bÞ
xy ¼ D u þ iqv ; r
ð27cÞ
e ¼ Dv þ iq u ;
ð28Þ
fD4 ½s2 ð1 þ s þ emÞ þ sð1 þ eÞ þ 2q2 D2 þ ss4 þ s3 þ s2 q2 ð1 þ s þ emÞ þ sq2 ð1 þ eÞ þ q4 g h ¼ 0:
ð29Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where D = (o/oy), and a ¼ b2 s2 þ q2 . Eliminating e between Eqs. (23) and (24), we obtain
The solution of Eq. (29), which is bounded at infinity, can be written as h ¼ Aðk 21 q2 s2 Þek1 y þ Bðk 22 q2 s2 Þek2 y ;
ð30Þ
178
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
where A and B are some parameters depending on q and s, and k1 and k2 are the roots with positive real parts of the characteristic equation k 4 ½s2 ð1 þ s þ emÞ þ sð1 þ eÞ þ 2q2 k 2 þ ss4 þ s3 þ s2 q2 ð1 þ s þ emÞ þ sq2 ð1 þ eÞ þ q4 ¼ 0;
ð31Þ
k1 and k2 are given by 1 2 k 21;2 ¼ q2 þ s ð1 þ s þ emÞ þ sð1 þ eÞ s 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ½sð1 þ s þ emÞ þ 1 þ e 4sð1 þ ssÞ :
ð32Þ
Substituting from Eq. (30) into Eq. (23), we obtain 1 e ¼ fAðk 21 q2 s ss2 Þðk 21 q2 s2 Þek1 y es þ Bðk 22 q2 s ss2 Þðk 22 q2 s2 Þek2 y g: ð33Þ and e satisfy Eq. (24), we obtain the compatibility Now, since both h conditions ðk 21 q2 s ss2 Þðk 21 q2 s2 Þ ¼ esð1 þ msÞðk 21 q2 Þ
ð34Þ
ðk 22 q2 s ss2 Þðk 22 q2 s2 Þ ¼ esð1 þ msÞðk 22 q2 Þ:
ð35Þ
and Using Eqs. (34) and (35), we can write Eq. (33) in the simplified form e ¼ ð1 þ msÞ½Aðk 21 q2 Þek1 y þ Bðk 22 q2 Þek2 y :
ð36Þ
Substituting from Eqs. (30) and (36) into Eq. (25), then solving the resulting equation, we obtain u ¼ iqð1 þ msÞ½Ceay þ Aek1 y þ Bek2 y ;
ð37Þ
where C is a parameter depending on q and s. Substituting from Eqs. (30) and (36) into Eq. (28) and solving the resulting equation, we get 2
q v ¼ ð1 þ msÞ Ceay þ k 1 Aek1 y þ k 2 Bek2 y : ð38Þ a Substituting from Eqs. (30) and (36)–(38) into Eq. (27), we obtain xx ¼ ð1 þ msÞ½ðb2 s2 2k 21 ÞAek1 y þ ðb2 s2 2k 22 ÞBek2 y 2q2 Ceay ; r ð39aÞ yy ¼ ð1 þ msÞ½ða2 þ q2 ÞAek1 y þ ða2 þ q2 ÞBek2 y þ 2q2 Ceay ; r
ð39bÞ
1 xy ¼ iqð1 þ msÞ 2k 1 Aek1 y þ 2k 2 Bek2 y þ ða2 þ q2 ÞCeay : r a
ð39cÞ
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
179
Taking the Laplace and Fourier Transform of Eqs. (14) and (15), respectively, we obtain the boundary conditions in the transformed domain as yy ðq; 0; sÞ ¼ f ðq; sÞ; r
ð40Þ
xy ðq; 0; sÞ ¼ 0; r
ð41Þ
ðq; sÞ: h ðq; 0; sÞ ¼ G
ð42Þ
Eqs. (30), (39a) and (39b) together with Eqs. (40)–(42) give
Aðk 21 q2 s2 Þ þ Bðk 22 q2 s2 Þ ¼ f ðq; sÞ;
ð43Þ
ðq; sÞ G ; ð1 þ msÞ
ð44Þ
ða2 þ q2 ÞA þ ða2 þ q2 ÞB þ 2q2 C ¼
2k 1 A þ 2k 2 B þ
ða2 þ q2 Þ C ¼ 0: a
ð45Þ
4. Inversion of the double transform We shall now outline the numerical inversion method used to find the solu tion in the physical domain. Let f ðq; y; sÞ be the Laplace–Fourier transform of a function f(x, y, t). First, we use the inversion formula of Fourier transform mentioned earlier to obtain a Laplace expression f ðx; y; sÞ of the form Z 1 1 eiqy f ðq; y; sÞ dq f ðx; y; sÞ ¼ pffiffiffiffiffiffi 2p 1 rffiffiffi Z 1 2 ¼ ðcosðqxÞf e ðq; y; sÞ þ sinðqxÞf o ðq; y; sÞÞ dq; p 0
where f e and f o denote the even and odd parts of f ðq; y; sÞ respectively. The complex inversion formula for Laplace transforms can be written as [16] Z dþi1 1 est f ðx; y; sÞ ds; f ðx; y; tÞ ¼ 2pi di1 where d is an arbitrary real number greater than all the real parts of the singularities of f ðq; y; sÞ. Taking s = d + ir, the above integral takes the form Z edt 1 itr f ðx; y; tÞ ¼ e f ðx; y; c þ irÞ dr: 2p 1 Expanding the function h(x, y, t) = exp(dt)f(x, y, t) in a Fourier series in the interval [0, 2L], we obtain the approximate formula [16]
180
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
f ðx; y; tÞ ¼ f1 ðx; y; tÞ þ ED ; where 1 X 1 f1 ðx; y; tÞ ¼ c0 ðx; y; tÞ þ ck ðx; y; tÞ for 0 6 t 6 2L 2 k¼1
ð46Þ
and ck ðx; y; tÞ ¼
edt Re½eikpt=L f ðx; y; d þ ikp=LÞ ; L
k ¼ 0; 1; 2; . . .
ð47Þ
ED, the discretization error, can be made arbitrarily small by choosing the constant d large enough [16]. Since the infinite series in Eq. (46) can only be summed up to a finite number N of terms, the approximate value of f(x, y, t) becomes N X 1 fN ðx; y; tÞ ¼ c0 ðx; y; tÞ þ ck ðx; y; tÞ 2 k¼1
for 0 6 t 6 2L:
ð48Þ
Using the above formula to evaluate f(x, y, 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 [16] 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(x, y, t): f ðx; y; tÞ ¼ f1 ðx; y; tÞ e2 dL f1 ðx; y; 2L þ tÞ þ E0D ; where the discretization error jE0D j jED j. Thus, the approximate value of f(x, y, t) becomes fNK ðx; y; tÞ ¼ fN ðx; y; tÞ e2 dL fN 0 ðx; y; 2L þ tÞ;
ð49Þ
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. (48). Let N be an odd natural number, and let Xm sm ðx; y; tÞ ¼ ck ðx; y; tÞ k¼1 be the sequence of partial sums of (48). We define the e-sequence by e0;m ¼ 0;
e1;m ¼ sm
and epþ1;m ¼ ep1;mþ1 þ 1=ðep;mþ1 ep;m Þ; p ¼ 1; 2; 3; . . .
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
181
It can be shown that [16] the sequence e1;1 ; e3;1 ; e5;1 ; . . . ; eN ;1 converges to f(q, y, 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. (49) together with the e-algorithm. The values of d and L are chosen according to the criteria outlined in [16].
5. Numerical results The functions f(x, t) and G(x, t) which represent the thermal and mechanical shock respectively are taken as f ðx; tÞ ¼ h0 H ða jxjÞ;
ð50Þ
Gðx; tÞ ¼ r0 H ða jxjÞ;
ð51Þ
where H is the Heaviside unit step function, a, r0 and h0 are constants. This means that the applied thermomechanical shock acts only on a band of width 2a centered around the x-axis on the surface of the half space and zero everywhere else. Then the applied thermomechanical shock in the transformed Laplace Fourier domains is given by the following two equations rffiffiffi 2 h0 sin qa ; ð52Þ f ðq; sÞ ¼ p qs rffiffiffi ðq; sÞ ¼ 2 r0 sin qa : G ð53Þ p qs The copper material was chosen for purposes of numerical evaluations. The constants of the problem are shown in Table 1. The computations were performed for three values of non-dimensional time, namely t = 0.04, t = 0.08 and t = 0.12. The numerical technique outlined above was used to obtain the temperature, displacement, and stress distributions. All
Table 1 Constants of the problem at = 1.78 (10)5 k = 7.76 (10)10 T0 = 293 a=1
k = 386 l = 3.86 (10)10 b2 = 4 h0 = 1
cE = 383.1 q = 8954 e = 0.0168 m = 0.02
g = 8886.73 c1 = 4.158 (10)3 s0 = 0.02 r0 = 1
182
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
1.20
θ
1.00
t =0.04 0.80
t =0.08 0.60
t =0.12 0.40 0.20 0.00 0.00
0.20
0.40
0.60
0.80
1.00
y
Temperature Distribution Fig. 1. Temperature distribution.
functions were evaluated inside the medium on the y-axis (x = 0) as functions of y. The temperature increment h is represented by the graph in Fig. 1. The displacement component v is shown in Fig. 2. The stress component ryy and rxx are shown in Figs. 3 and 4. The displacement is continuous for all values of y. The temperature has one discontinuity at y = 0.28013 for t = 0.04, at y = 0.56149 for t = 0.08 and at y = 0.8422 for t = 0.12. Both normal and tangential stress distributions has two discontinuities at y = 0.0317 and y = 0.28013 for t = 0.04, at y = 0.0724 and y = 0.56149 for t = 0.08 and at y = 0.1106 and y = 0.8422 for t = 0.12.
0.04
v
0.02 0.00
t = 0.04
-0.02
t = 0.08
-0.04
t = 0.12
-0.06 -0.08 0
0.2
0.4
0.6
0.8
y
Normal Displacement Distribution Fig. 2. Normal displacement distribution.
1
1.2
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
183
σ yy
0.0
-1.0
t = 0.04 -2.0
t = 0.08 t = 0.12
-3.0
-4.0 0
0.2
0.4
0.6
0.8
1
1.2
y
Normal Stress Distribution Fig. 3. Normal stress distribution.
σ xx
0.0
-1.0
t = 0.04 -2.0
t = 0.08 t =0.12
-3.0
-4.0 0
0.2
0.4
0.6
0.8
1
1.2
y
Tangential Stress Distribution Fig. 4. Tangential stress distribution.
References [1] V. Danilovskaya, Thermal stresses in an elastic half-space due to a sudden heating of its boundary, Prikl. Mat. Mekh. 14 (1950) 316–321 (In Russian). [2] M. Biot, Thermoelasticity and irreversible thermo-dynamics, J. Appl. Phys. 27 (1956) 240–253. [3] H. Lord, Y. Shulman, A generalized dynamical theory of thermoelasticity, J. Mech. Phys. Solids 15 (1967) 299–309. [4] R. Dhaliwal, H. Sherief, Generalized thermoelasticity for anisotropic media, Quart. Appl. Math. 33 (1980) 1–8. [5] J. Ignaczak, Uniqueness in generalized thermoelasticity, J. Thermal Stresses 2 (1979) 171–175. [6] J. Ignaczak, A note on uniqueness in thermoelasticity with one relaxation time, J. Thermal Stresses 5 (1982) 257–263.
184
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
[7] H. Sherief, R. Dhaliwal, A uniqueness theorem and a variational principal for generalized thermoelasticity, J. Thermal Stresses 3 (1980) 223–230. [8] H. Sherief, On uniqueness and stability in generalized thermoelasticity, Quart. Appl. Math. 45 (1987) 773–778. [9] A.E. Green, K.A. Lindsay, Thermoelasticity, J. Elast. 2 (1972) 1–7. [10] A.E. Green, A note in linear, Mathematika 19 (1972) 69–75. [11] H. Sherief, Fundamental solution of the generalized thermoelastic problem for short times, J. Thermal Stresses 9 (1986) 151–164. [12] H. Sherief, K. Helmy, A two dimensional generalized thermoelasticity problem for a halfspace, J. Thermal Stresses 22 (1999) 897–910. [13] H. Sherief, N.M. El-Maghraby, An internal penny shaped crack in an infinite thermoelastic soiled, J. Thermal Stresses 26 (2003) 333–352. [14] N.M. El-Maghraby, Two dimensional problem in generalized thermoelasticity with heat sources, J. Thermal Stresses 27 (2004) 227–239. [15] N.M. El-Maghraby, Hamdy M. Yossef, State space approach to generalized thermoelastic problem with thermomechanical shock, Appl. Math. Comput. 156 (2004) 577–586. [16] G. Honig, U. Hirdes, A method for the numerical inversion of the Laplace transform, J. Comp. Appl. Math. 10 (1984) 113–132.
A two-dimensional thermoelasticity problem for thermomechanical shock with two relaxation times Nasser M. El-Maghraby *, Hamdy M. Yossef Department of Basic and Applied Science, Arab Academy for Science and Technology, P.O. Box 1029, Alexandria, Egypt
Abstract In this work a two-dimensional problem of thermoelasticity with two relaxation times is introduced. Laplace and Fourier transform techniques are used. The resulting formulation is applied to a thermomechanical shock half-space problem. The solution in the transformed domain obtained by using a direct approach. Numerical inversion of both transforms is carried out to obtain the temperature, stress and displacement distributions in the physical domain. Numerical results are represented graphically. Ó 2004 Elsevier Inc. All rights reserved. Keywords: Elasticity; Thermoelasticity; Two relaxation times; Mechanical shock
1. Introduction Danilovskaya [1] was the first to solve an actual problem in the theory of elasticity with uniform heat. The problem was for a half-space subjected to a *
Corresponding author. E-mail address: [email protected] (N.M. El-Maghraby).
0096-3003/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.10.077
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
173
thermal shock in the context of what became known as the theory of uncoupled thermoelasticity. In this theory the temperature is governed by a parabolic partial differential equation which does not contain any elastic terms. It was not much latter that many attempts were made to remedy the shortcomings of this theory. Biot [2] formulated the theory of coupled thermoelasticity to eliminate the paradox inherent in the classical uncoupled theory that elastic changes have no effect on the temperature. The heat equations for both theories, however, are of the diffusion type predicting infinite speeds of propagation for heat waves contrary to physical observations. Lord and Shulman [3] introduced the theory of generalized thermoelasticity with one relaxation time for the special case of an isotropic body. This theory was extended by Dhaliwal and Sherief [4] to include the anisotropic case. In this theory a modified law of heat conduction including both the heat flux and its time derivative replaces the conventional FourierÕs law. The heat equation associated with this theory is hyperbolic and hence eliminates the paradox of infinite speeds of propagation inherent in both the uncoupled and coupled theories of thermoelasticity. Uniqueness of solution for this theory was proved under different conditions by Ignaczak [5,6], Sherief and Dhaliwal [7], Dhaliwal and Sherief [4] and Sherief [8]. Green and Lindsay [9] obtained the theory of generalized thermoelasticity with two relaxation times. The uniqueness of solution for this theory was proved by Green [10]. The fundamental solution was obtained by Sherief [11]. Sherief and Helmy have solved a two-dimensional problem in [12]. Sherief and El-Maghraby have solved a crack problem in one relaxation time [13], ElMaghraby [14] has solved a two-dimensional problem for a half space subjected to heat sources. The authors have solved a one-dimensional problem with thermomechanical shock [15]. In this work a two-dimensional problem of thermoelasticity with two relaxation times is introduced. Laplace and Fourier transform techniques are used. The resulting formulation is applied to a thermomechanical shock half-space problem. The solution in the transformed domain obtained by using a direct approach. Numerical inversion of both transforms is carried out to obtain the temperature, stress and displacement distributions in the physical domain. Numerical results are represented graphically.
2. Formulation of the problem We shall consider a homogeneous isotropic thermoelastic solid occupying the half-space y P 0. The y-axis is taken perpendicular to the bounding plane pointing inwards. We shall also assume that the initial state of the medium is quiescent. The surface of this medium is subjected to a thermomechanical shock. The displacement vector thus, has components.
174
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
u ¼ ðu; v; 0Þ and the cubical dilatation e is given by e ¼ div u ¼
ou ov þ : ox oy
ð1Þ
The governing equations for thermoelasticity with two relaxation times consist of: (1) The equation of motion in vector form o2 u oT q 2 ¼ ðk þ lÞ grad e þ lr2 u c grad T þ m ; ð2Þ ot ot which has two Cartesian components o2 u oe o2 u o2 u oT o2 T þm ; q 2 ¼ ðk þ lÞ þ l 2 þ l 2 c ot ox ox oy ox ox ot
ð3Þ
o2 v oe o2 v o2 v oT o2 T þ l þ m ¼ ðk þ lÞ þ l c ; ot2 oy ox2 oy 2 oy oy ot
ð4Þ
q
where k, l are Lame´Õs constants, q is the density, t is the time variable, T is the absolute temperature of the medium, m is a constant with dimensions of time, called a relaxation time, and c is a material constant given by c = (3k + 2l)at, at being the coefficients of linear thermal expansion. (2) The equation of heat conduction oT o2 T oe þ s 2 þ cT 0 ; ð5Þ kr2 T ¼ qC E ot ot ot where k is the thermal conductivity of the medium, T0 is a reference temperature, CE is the specific heat at constant strain and s is another relaxation time. Applying the div operator to both sides of Eq. (2), we obtain o2 e oT q 2 ¼ ðk þ 2lÞr2 e cr2 T þ m ; ð6Þ ot ot where $2 is LaplaceÕs operator. The above equations are supplemented with the constitutive equations ov oT rxx ¼ ðk þ 2lÞe 2l c T T 0 þ m ; ð7aÞ oy ot ou oT ryy ¼ ðk þ 2lÞe 2l c T T 0 þ m ; ð7bÞ ox ot
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
oT rzz ¼ ke c T T 0 þ m ; ot ou ov þ rxx ¼ l ; oy ox
175
ð7cÞ
ð7dÞ
rxz ¼ ryx ¼ 0:
ð7eÞ
We shall use the following non-dimensional variables x0 ¼ c1 gx;
y 0 ¼ c1 gy;
t0 ¼ c21 gt;
u0 ¼ c1 gu;
v0 ¼ c1 gv;
cðT T 0 Þ rij ; rij ¼ ; h¼ k þ 2l l pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where g = (qcE/k) and c1 ¼ ðk þ 2l=qÞ is the speed of propagation of isothermal elastic waves. Eq. (1) retains its form, while Eqs. (3)–(6) take the following form (dropping the primes for convenience). 2 oe o2 h 2o u 2 2 oh 2 þm b 2 ¼ ðb 1Þ þ r u b ; ð8Þ ot ox ox ox ot s0 ¼ c21 gs;
b2
m0 ¼ c21 gm;
o2 v oe o2 h 2 2 oh 2 þ r þ m ¼ ðb 1Þ v b ; ot2 oy oy ox ot
oh o2 h oe þs 2 þe ; ot ot ot o2 e oh b 2 2 ¼ r 2 e b2 r 2 h þ m ; ot ot r2 h ¼
ð9Þ
ð10Þ ð11Þ
where b2 ¼
k þ 2l l
and
e¼
T 0 c2 : qcE ðk þ 2lÞ
The constitutive equation (7) in non-dimensional form can be written as ov oh rxx ¼ b2 e 2 b2 h þ m ; ð12aÞ oy ot ou oh ryy ¼ b2 e 2 b2 h þ m ; ð12bÞ ox ot oh rzz ¼ ðb2 2Þe b2 h þ m ; ð12cÞ ot
176
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
rxy ¼
ou ov þ : oy ox
ð12dÞ
The above equations are solved subject to the initial conditions h ¼ u ¼ v ¼ 0;
oh ou ov ¼ ¼ ¼ 0 at t ¼ 0 ot ot ot
ð13Þ
and the boundary conditions rxy ¼ 0;
ryy ¼ f ðx; tÞ
h ¼ Gðx; tÞ
at y ¼ 0;
at y ¼ 0;
ð14Þ ð15Þ
where f(x, t) and G(x, t) are known functions. Using FourierÕs law of heat conduction in non-dimensional form (which is valid for Green theory), namely qn ¼
oh ; on
ð16Þ
where qn denotes the normal component of the heat flux vector.
3. Solution in the transformed domain Taking the Laplace transform with parameter s (denoted by a bar) of both sides of Eqs. (8)–(12), we obtain the following set of equations ðr2 s ss2 Þ h ¼ ese;
ð17Þ
ðr2 sÞe ¼ ð1 þ msÞr2 h;
ð18Þ
b 2 s2 u ¼ ðb2 1Þ
oe oh þ r2 u b2 ð1 þ msÞ ; ox ox
ð19Þ
b2 s2v ¼ ðb2 1Þ
oe o h þ r2v b2 ð1 þ msÞ ; oy oy
ð20Þ
xx ¼ b2e 2 r
ov b2 ð1 þ msÞ h; oy
ð21aÞ
yy ¼ b2e 2 r
o u b2 ð1 þ msÞ h; ox
ð21bÞ
xy ¼ r
o u ov þ oy ox
ð21cÞ
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
177
and the Laplace transform of the cubical dilation (1) becomes o u ov e ¼ þ : ð22Þ oy ox We have assumed that the Laplace transforms of the temperature, stress components, displacements, and the applied thermomechanical shock exist. We now introduce the exponential Fourier transform (denoted by an asterisk) with respect to the space variable x, defined by Z 1 1 f ðq; y; sÞ ¼ I½f ðx; y; sÞ ¼ pffiffiffiffiffiffi eiqx f ðx; y; sÞ dx 2p 1 with corresponding inversion formula Z 1 1 eiqx f ðq; y; sÞ dq f ðx; y; sÞ ¼ I1 ½f ðq; y; sÞ ¼ pffiffiffiffiffiffi 2p 1
where i ¼
pffiffiffiffiffiffiffi 1:
We assume that all the relevant functions (temperature, stress, . . ., etc.) are sufficient smooth on the real line such that the Fourier transforms of these functions exists. Taking the Fourier transform of both sides of Eqs. (17)–(22), we obtain the following set of equations ðD2 q2 s ss2 Þ h ¼ ese ; ð23Þ ðD2 q2 s2 Þe ¼ ð1 þ msÞðD2 q2 Þ h;
ð24Þ
u ¼ iq½b2 ð1 þ msÞ h ðb2 1Þe ; ðD2 a2 Þ
ð25Þ
h ðb2 1ÞDe ; ðD2 a2 Þv ¼ b2 ð1 þ msÞD
ð26Þ
xx ¼ b2e 2Dv b2 ð1 þ msÞ h; r
ð27aÞ
yy ¼ b2e 2iq r u b2 ð1 þ msÞ h;
ð27bÞ
xy ¼ D u þ iqv ; r
ð27cÞ
e ¼ Dv þ iq u ;
ð28Þ
fD4 ½s2 ð1 þ s þ emÞ þ sð1 þ eÞ þ 2q2 D2 þ ss4 þ s3 þ s2 q2 ð1 þ s þ emÞ þ sq2 ð1 þ eÞ þ q4 g h ¼ 0:
ð29Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where D = (o/oy), and a ¼ b2 s2 þ q2 . Eliminating e between Eqs. (23) and (24), we obtain
The solution of Eq. (29), which is bounded at infinity, can be written as h ¼ Aðk 21 q2 s2 Þek1 y þ Bðk 22 q2 s2 Þek2 y ;
ð30Þ
178
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
where A and B are some parameters depending on q and s, and k1 and k2 are the roots with positive real parts of the characteristic equation k 4 ½s2 ð1 þ s þ emÞ þ sð1 þ eÞ þ 2q2 k 2 þ ss4 þ s3 þ s2 q2 ð1 þ s þ emÞ þ sq2 ð1 þ eÞ þ q4 ¼ 0;
ð31Þ
k1 and k2 are given by 1 2 k 21;2 ¼ q2 þ s ð1 þ s þ emÞ þ sð1 þ eÞ s 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ½sð1 þ s þ emÞ þ 1 þ e 4sð1 þ ssÞ :
ð32Þ
Substituting from Eq. (30) into Eq. (23), we obtain 1 e ¼ fAðk 21 q2 s ss2 Þðk 21 q2 s2 Þek1 y es þ Bðk 22 q2 s ss2 Þðk 22 q2 s2 Þek2 y g: ð33Þ and e satisfy Eq. (24), we obtain the compatibility Now, since both h conditions ðk 21 q2 s ss2 Þðk 21 q2 s2 Þ ¼ esð1 þ msÞðk 21 q2 Þ
ð34Þ
ðk 22 q2 s ss2 Þðk 22 q2 s2 Þ ¼ esð1 þ msÞðk 22 q2 Þ:
ð35Þ
and Using Eqs. (34) and (35), we can write Eq. (33) in the simplified form e ¼ ð1 þ msÞ½Aðk 21 q2 Þek1 y þ Bðk 22 q2 Þek2 y :
ð36Þ
Substituting from Eqs. (30) and (36) into Eq. (25), then solving the resulting equation, we obtain u ¼ iqð1 þ msÞ½Ceay þ Aek1 y þ Bek2 y ;
ð37Þ
where C is a parameter depending on q and s. Substituting from Eqs. (30) and (36) into Eq. (28) and solving the resulting equation, we get 2
q v ¼ ð1 þ msÞ Ceay þ k 1 Aek1 y þ k 2 Bek2 y : ð38Þ a Substituting from Eqs. (30) and (36)–(38) into Eq. (27), we obtain xx ¼ ð1 þ msÞ½ðb2 s2 2k 21 ÞAek1 y þ ðb2 s2 2k 22 ÞBek2 y 2q2 Ceay ; r ð39aÞ yy ¼ ð1 þ msÞ½ða2 þ q2 ÞAek1 y þ ða2 þ q2 ÞBek2 y þ 2q2 Ceay ; r
ð39bÞ
1 xy ¼ iqð1 þ msÞ 2k 1 Aek1 y þ 2k 2 Bek2 y þ ða2 þ q2 ÞCeay : r a
ð39cÞ
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
179
Taking the Laplace and Fourier Transform of Eqs. (14) and (15), respectively, we obtain the boundary conditions in the transformed domain as yy ðq; 0; sÞ ¼ f ðq; sÞ; r
ð40Þ
xy ðq; 0; sÞ ¼ 0; r
ð41Þ
ðq; sÞ: h ðq; 0; sÞ ¼ G
ð42Þ
Eqs. (30), (39a) and (39b) together with Eqs. (40)–(42) give
Aðk 21 q2 s2 Þ þ Bðk 22 q2 s2 Þ ¼ f ðq; sÞ;
ð43Þ
ðq; sÞ G ; ð1 þ msÞ
ð44Þ
ða2 þ q2 ÞA þ ða2 þ q2 ÞB þ 2q2 C ¼
2k 1 A þ 2k 2 B þ
ða2 þ q2 Þ C ¼ 0: a
ð45Þ
4. Inversion of the double transform We shall now outline the numerical inversion method used to find the solu tion in the physical domain. Let f ðq; y; sÞ be the Laplace–Fourier transform of a function f(x, y, t). First, we use the inversion formula of Fourier transform mentioned earlier to obtain a Laplace expression f ðx; y; sÞ of the form Z 1 1 eiqy f ðq; y; sÞ dq f ðx; y; sÞ ¼ pffiffiffiffiffiffi 2p 1 rffiffiffi Z 1 2 ¼ ðcosðqxÞf e ðq; y; sÞ þ sinðqxÞf o ðq; y; sÞÞ dq; p 0
where f e and f o denote the even and odd parts of f ðq; y; sÞ respectively. The complex inversion formula for Laplace transforms can be written as [16] Z dþi1 1 est f ðx; y; sÞ ds; f ðx; y; tÞ ¼ 2pi di1 where d is an arbitrary real number greater than all the real parts of the singularities of f ðq; y; sÞ. Taking s = d + ir, the above integral takes the form Z edt 1 itr f ðx; y; tÞ ¼ e f ðx; y; c þ irÞ dr: 2p 1 Expanding the function h(x, y, t) = exp(dt)f(x, y, t) in a Fourier series in the interval [0, 2L], we obtain the approximate formula [16]
180
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
f ðx; y; tÞ ¼ f1 ðx; y; tÞ þ ED ; where 1 X 1 f1 ðx; y; tÞ ¼ c0 ðx; y; tÞ þ ck ðx; y; tÞ for 0 6 t 6 2L 2 k¼1
ð46Þ
and ck ðx; y; tÞ ¼
edt Re½eikpt=L f ðx; y; d þ ikp=LÞ ; L
k ¼ 0; 1; 2; . . .
ð47Þ
ED, the discretization error, can be made arbitrarily small by choosing the constant d large enough [16]. Since the infinite series in Eq. (46) can only be summed up to a finite number N of terms, the approximate value of f(x, y, t) becomes N X 1 fN ðx; y; tÞ ¼ c0 ðx; y; tÞ þ ck ðx; y; tÞ 2 k¼1
for 0 6 t 6 2L:
ð48Þ
Using the above formula to evaluate f(x, y, 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 [16] 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(x, y, t): f ðx; y; tÞ ¼ f1 ðx; y; tÞ e2 dL f1 ðx; y; 2L þ tÞ þ E0D ; where the discretization error jE0D j jED j. Thus, the approximate value of f(x, y, t) becomes fNK ðx; y; tÞ ¼ fN ðx; y; tÞ e2 dL fN 0 ðx; y; 2L þ tÞ;
ð49Þ
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. (48). Let N be an odd natural number, and let Xm sm ðx; y; tÞ ¼ ck ðx; y; tÞ k¼1 be the sequence of partial sums of (48). We define the e-sequence by e0;m ¼ 0;
e1;m ¼ sm
and epþ1;m ¼ ep1;mþ1 þ 1=ðep;mþ1 ep;m Þ; p ¼ 1; 2; 3; . . .
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
181
It can be shown that [16] the sequence e1;1 ; e3;1 ; e5;1 ; . . . ; eN ;1 converges to f(q, y, 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. (49) together with the e-algorithm. The values of d and L are chosen according to the criteria outlined in [16].
5. Numerical results The functions f(x, t) and G(x, t) which represent the thermal and mechanical shock respectively are taken as f ðx; tÞ ¼ h0 H ða jxjÞ;
ð50Þ
Gðx; tÞ ¼ r0 H ða jxjÞ;
ð51Þ
where H is the Heaviside unit step function, a, r0 and h0 are constants. This means that the applied thermomechanical shock acts only on a band of width 2a centered around the x-axis on the surface of the half space and zero everywhere else. Then the applied thermomechanical shock in the transformed Laplace Fourier domains is given by the following two equations rffiffiffi 2 h0 sin qa ; ð52Þ f ðq; sÞ ¼ p qs rffiffiffi ðq; sÞ ¼ 2 r0 sin qa : G ð53Þ p qs The copper material was chosen for purposes of numerical evaluations. The constants of the problem are shown in Table 1. The computations were performed for three values of non-dimensional time, namely t = 0.04, t = 0.08 and t = 0.12. The numerical technique outlined above was used to obtain the temperature, displacement, and stress distributions. All
Table 1 Constants of the problem at = 1.78 (10)5 k = 7.76 (10)10 T0 = 293 a=1
k = 386 l = 3.86 (10)10 b2 = 4 h0 = 1
cE = 383.1 q = 8954 e = 0.0168 m = 0.02
g = 8886.73 c1 = 4.158 (10)3 s0 = 0.02 r0 = 1
182
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
1.20
θ
1.00
t =0.04 0.80
t =0.08 0.60
t =0.12 0.40 0.20 0.00 0.00
0.20
0.40
0.60
0.80
1.00
y
Temperature Distribution Fig. 1. Temperature distribution.
functions were evaluated inside the medium on the y-axis (x = 0) as functions of y. The temperature increment h is represented by the graph in Fig. 1. The displacement component v is shown in Fig. 2. The stress component ryy and rxx are shown in Figs. 3 and 4. The displacement is continuous for all values of y. The temperature has one discontinuity at y = 0.28013 for t = 0.04, at y = 0.56149 for t = 0.08 and at y = 0.8422 for t = 0.12. Both normal and tangential stress distributions has two discontinuities at y = 0.0317 and y = 0.28013 for t = 0.04, at y = 0.0724 and y = 0.56149 for t = 0.08 and at y = 0.1106 and y = 0.8422 for t = 0.12.
0.04
v
0.02 0.00
t = 0.04
-0.02
t = 0.08
-0.04
t = 0.12
-0.06 -0.08 0
0.2
0.4
0.6
0.8
y
Normal Displacement Distribution Fig. 2. Normal displacement distribution.
1
1.2
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
183
σ yy
0.0
-1.0
t = 0.04 -2.0
t = 0.08 t = 0.12
-3.0
-4.0 0
0.2
0.4
0.6
0.8
1
1.2
y
Normal Stress Distribution Fig. 3. Normal stress distribution.
σ xx
0.0
-1.0
t = 0.04 -2.0
t = 0.08 t =0.12
-3.0
-4.0 0
0.2
0.4
0.6
0.8
1
1.2
y
Tangential Stress Distribution Fig. 4. Tangential stress distribution.
References [1] V. Danilovskaya, Thermal stresses in an elastic half-space due to a sudden heating of its boundary, Prikl. Mat. Mekh. 14 (1950) 316–321 (In Russian). [2] M. Biot, Thermoelasticity and irreversible thermo-dynamics, J. Appl. Phys. 27 (1956) 240–253. [3] H. Lord, Y. Shulman, A generalized dynamical theory of thermoelasticity, J. Mech. Phys. Solids 15 (1967) 299–309. [4] R. Dhaliwal, H. Sherief, Generalized thermoelasticity for anisotropic media, Quart. Appl. Math. 33 (1980) 1–8. [5] J. Ignaczak, Uniqueness in generalized thermoelasticity, J. Thermal Stresses 2 (1979) 171–175. [6] J. Ignaczak, A note on uniqueness in thermoelasticity with one relaxation time, J. Thermal Stresses 5 (1982) 257–263.
184
N.M. El-Maghraby, H.M. Yossef / Appl. Math. Comput. 170 (2005) 172–184
[7] H. Sherief, R. Dhaliwal, A uniqueness theorem and a variational principal for generalized thermoelasticity, J. Thermal Stresses 3 (1980) 223–230. [8] H. Sherief, On uniqueness and stability in generalized thermoelasticity, Quart. Appl. Math. 45 (1987) 773–778. [9] A.E. Green, K.A. Lindsay, Thermoelasticity, J. Elast. 2 (1972) 1–7. [10] A.E. Green, A note in linear, Mathematika 19 (1972) 69–75. [11] H. Sherief, Fundamental solution of the generalized thermoelastic problem for short times, J. Thermal Stresses 9 (1986) 151–164. [12] H. Sherief, K. Helmy, A two dimensional generalized thermoelasticity problem for a halfspace, J. Thermal Stresses 22 (1999) 897–910. [13] H. Sherief, N.M. El-Maghraby, An internal penny shaped crack in an infinite thermoelastic soiled, J. Thermal Stresses 26 (2003) 333–352. [14] N.M. El-Maghraby, Two dimensional problem in generalized thermoelasticity with heat sources, J. Thermal Stresses 27 (2004) 227–239. [15] N.M. El-Maghraby, Hamdy M. Yossef, State space approach to generalized thermoelastic problem with thermomechanical shock, Appl. Math. Comput. 156 (2004) 577–586. [16] G. Honig, U. Hirdes, A method for the numerical inversion of the Laplace transform, J. Comp. Appl. Math. 10 (1984) 113–132.