10_xxx_kar.dvi EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 4, No. 3, 2011, 304-321 ISSN 1307-5543 – www.ejpam.com Analysis of Thermoelastic Response in a Fiber Reinforced Thin Annular Disc With Three-Phase-Lag Effect Avijit Kar1,∗, M. Kanoria2 1 River Research Institute, Haringhata Central Laboratory, Mohanpur, Nadia-741246, India 2 Department of Applied Mathematics, University of Calcutta, 92, A.P.C. Road, Kolkata - 700 009, India Abstract. This paper concerns the investigation of the stresses, displacement and temperature due to the axisymmetric thermal shock loading on the inner boundary in a transversely isotropic thin annular disc. The analysis encompasses thermo-elasticity without energy dissipation theory (TEWOED (GN- II)) and three-phase-lag theory of the generalized thermo-elasticity to account for the finite velocity of the temperature. The Laplace transform method is used to transform the coupled equations into the space domain, where two different methods, eigen-value approach and the Galerkin finite element technique are employed to solve the resulting equations in the transformed domain. The inverse of the transformed solution is carried out by applying a method of Bellman et al. Stresses, displacement and temperature distributions have been computed numerically and presented graphically in a number of figures. A comparison of the results for different theories (GN-II and Three-phase-lag model) and for two different methods are presented. 2000 Mathematics Subject Classifications: 74F Key Words and Phrases: Generalized thermo-elasticity, energy dissipation, Laplace transform, eigen- value approach, Galerkin finite element technique. 1. Introduction Generalized thermo-elasticity theories have been developed with the objective of removing the paradox of infinite speed of thermal signals inherent in the conventional coupled dynami- cal theory of thermo-elasticity in which parabolic type heat conduction equation is considered, contradict physical facts. During the last three decades, generalized theories involving finite speed of heat transportation (hyperbolic heat transport equation) in elastic solids have been developed to remove this paradox. The first generalization is proposed by Lord and Shulman [23] and is known as extended thermo-elasticity theory (ETE), which involves one thermal re- laxation time parameter (single-phase-lag model). The second generalization to the coupled ∗Corresponding author. Email addresses: akmathemati s�yahoo. om (A. Kar), k_mri�yahoo. om (M. Kanoria) http://www.ejpam.com 304 c© 2011 EJPAM All rights reserved. A. Kar, M. Kanoria / Eur. J. Pure Appl. Math, 4 (2011), 304-321 305 thermo-elasticity theory is developed by Green and Lindsay [12], which involving two relax- ation times is known as temperature rate dependent thermo-elasticity (TRDTE). Experimental studies indicate that the relaxation times can be of relevance in the cases involving a rapidly propagating crack tip, shock waves propagation, laser technique etc. Because of the experi- mental evidence in support of finiteness of heat propagation speed, the generalized thermo- elasticity theories are considered to be more realistic than the conventional theory in dealing with practical problems involving very large heat fluxes at short intervals like those occurring in laser units and energy channels. The third generalization is known as low-temperature thermo-elasticity introduced by Hetnarski and Ignaczak [16] called H-I theory. Most engi- neering materials such as metals possess a relatively high rate of thermal damping and thus are not suitable for use in experiments concerning second sound propagation. But, given the state of recent advances in material science, it may be possible in the foreseeable future to identify (or even manufacture for laboratory purposes) an idealized material for the pur- pose of studying the propagation of thermal waves at finite speed. The fourth generalization is concerned with the thermo-elasticity without energy dissipation (TEWOED) and thermo- elasticity with energy dissipation (TEWED) introduced by Green and Naghdi [13, 14, 15] and provides sufficient basic modifications in the constitutive equations that permit treatment of a much wider class of heat flow problems, labelled as types I, II, III. The natures of these three types of constitutive equations are such that when the respective theories are linearized, type- I is the same as the classical heat equation (based on FourierŠs law) whereas the linearized versions of type-II and type-III theories permit propagation of thermal waves at finite speed. The entropy flux vector in type II and III ( i.e. thermo-elasticity without energy dissipation (TEWOED) and thermo-elasticity with energy dissipation (TEWED)) models are determined in terms of the potential that also determines stresses. When Fourier conductivity is dominant the temperature equation reduces to classical Fourier law of heat conduction and when the effect of conductivity is negligible the equation has undamped thermal wave solutions with- out energy dissipation. Applying the above theories of generalized thermo-elasticity, several problems have been solved by Bagri and Eslami [1, 3, 2], Kar and Kanoria [17, 18, 19], Das et al. [8, 21], Roychoudhuri and Datta [29], Roychoudhury and Bandyopadhyay [28], Chan- drasekharaiah [5, 6], Wang et al. [32], Ghosh and Kanoria [10, 11], Mallik and Kanoria [24], Ootao et al. [25], Shao et al. [30], Wang and Mai [33] etc. The fifth generalization to the thermo-elasticity theory is known as the dual-phase-lag thermo- elasticity developed by Tzou [31] and Chandrasekhariah [7]. Tzou considered microstruc- tural effects into the delayed response in time in the macroscopic formulation by taking into account that increase of the lattice temperature is delayed due to photon-electron interactions on the macroscopic level. Tzou [31] introduced two-phase-lags to both the heat flux vector and the temperature gradient. According to this model, classical Fourier’s law −→q = −K −→∇ T has been replaced by −→q (P, t +τq) = −K −→∇ T (P, t +τT ), where the temperature gradient −→∇ T at a point P of the material at time t +τT corresponds to the heat flux vector −→q at the same point at time t + τq. Here K is the thermal conductivity of the material. The delay time τT is interpreted as that caused by the microstructural interactions and is called the phase-lag of the temperature gradient. The other delay time τq is interpreted as the relaxation time due A. Kar, M. Kanoria / Eur. J. Pure Appl. Math, 4 (2011), 304-321 306 to the fast transient effects of thermal inertia and is called the phase-lag of the heat flux. For τq = τT = 0, the Fourier’s law in two-phase-lag model is identical with classical Fourier’s law. If τq = τ and τT = 0, Tzou [31] refers to the model as single-phase-lag model. Roychoudhuri [26] studied one-dimensional thermo-elastic wave propagation in an elastic half-space in the context of dual-phase-lag model. The sixth generalization is known as three-phase-lag thermo-elasticity which is due to Roy- choudhuri [27]. According to this model−→q (P, t+τq) = −[K−→∇ T (P, t+τT )+K∗ −→∇ν(P, t+τν)], where −→∇ν (ν̇ = T) is the thermal displacement gradient and K∗ is the additional material constant and τν is the phase lag for thermal displacement gradient. To study some prac- tical relevant problems and have found that in heat transfer problems involving very short time intervals and in the problems of very high heat fluxes, the hyperbolic equation gives significantly different results than the parabolic equation. According to this phenomenon the lagging behavior in the heat conduction in solid should not be ignored particularly when the elapsed times during a transient process are very small, say about 10−7 second or the heat flux is very much high. Three-phase-lag model is very useful in the problems of nuclear boiling, exothermic catalytic reactions, phonon-electron interactions, phonon-scattering etc., where the delay time τq captures the thermal wave behavior (a small scale response in time), the phase-lag τT captures the effect of phonon-electron interactions (a microscopic response in space), the other delay time τν is effective since, in the three-phase-lag model, the thermal displacement gradient is considered as a constitutive variable whereas in the conventional thermo-elasticity theory temperature gradient is considered as a constitutive variable. Re- cently, Kar and Kanoria [20] studied thermo-visco-elastic problem of a spherical shell in the context of three-phase-lag model. However, over the last few decades anisotropic materials have been increasingly used. There are materials which have natural anisotropy such as zinc, magnesium, sapphire, wood, some rocks and crystals, and also there are artificially manufactured materials such as fibre-reinforced composite materials, which exhibit anisotropic character. The advantage of composite mate- rials over the traditional materials lies on their valuable strength, elastic and other properties [22]. A reinforced material may be regarded to some order of approximation, as homoge- neous and anisotropic elastic medium having a certain kind of elastic symmetry depending on the symmetry of reinforcement. Some glass fibre reinforced plastics may be regarded as trans- versely isotropic. Thus problems of solid mechanics should not be restricted to the isotropic medium only. Increasing use of anisotropic media demands that the study of elastic problems should be extended to anisotropic medium also. To the authors’ knowledge, under three-phase-lag effect no solution of transversely isotropic materials has been reported. With this motivation in mind the present analysis is to study the thermo-elastic stresses, displacement and temperature distribution in a transversely isotropic thin annular disc due to the axisymmetric thermal shock loading on the inner boundary of the disc in the context of TEWOED [14] of type-II and three-phase-lag model of generalized thermo-elasticity. The Laplace transformation method is used to transform the equations from the time domain to the Laplace domain. Two different methods, (A) eigen-value approach and (B) the Galerkin finite element technique are employed to solve the resulting equations in the A. Kar, M. Kanoria / Eur. J. Pure Appl. Math, 4 (2011), 304-321 307 transformed domain. In the first method (eigen value approach) the fundamental equations have been expressed in the form of vector-matrix differential equation which is then solved by eigen value approach [21] and in the second method (Galerkin finite element method) the radius of the disc is discretized into a finite elements along the radial direction where the Galerkin finite element method is employed to derive the force and stiffness matrices of the base element. Then the system of equations is solved to find the nodal values of the stresses, temperature and displacement. Finally, inversion of the Laplace transform is done following Bellman et al. [4]. The results obtained theoretically have been computed numerically and are presented graphically for Sapphire material. A complete and comprehensive analysis and comparison of the results for different theories (GN-II and Three-phase-lag model) as well as in two different methods are also presented. 2. Basic Equations and Constitutive Eelations We consider a homogeneous transversely isotropic thermo-elastic thin annular disc of in- ner radius a and outer radius b having initially undisturbed state at an uniform temperature T0, under axisymmetric thermal shock load applied into its inner boundary. We use plane polar co-ordinates (r,θ) with the centre of the hole as the origin. In the present problem (due to central symmetry) the displacement and temperature are as- sumed to be functions of r and time t only. The stress-strain-temperature relations in the generalized theory are [22] τr r = C11 ∂ u ∂ r + C12 u r −β11T, (1) τθθ = C12 ∂ u ∂ r + C11 u r − β22T (2) and the generalized heat conduction equation in the three-phase-lag model [27] is K∗∇2T +τ∗ν∇2 Ṫ + KτT∇2 T̈ = � 1+τq ∂ ∂ t + 1 2 τ2 q ∂ 2 ∂ t2 � × � ρCe T̈ + T0 ∂ 2 ∂ t2 � β11 ∂ u ∂ r + β22 u r �� , (3) where τi j (i, j = r,θ) is the stress tensor, T is the temperature increase over the reference temperature T0, C11 and C12 are elastic constants, β11 and β22 are thermal modulii, K is the coefficient of thermal conductivity along radial direction, K∗ is the additional material constant along radial direction, ρ is the mass density, Ce is the specific heat of the solid at constant strain, τT and τq are the the phase-lag of temperature gradient and the phase-lag of heat flux respectively. Also τ∗ν = K+τνK ∗, where τν is the phase-lag of thermal displacement gradient. Equations (1) - (3), when K = τT = τq = τν = 0, reduce to the equations of A. Kar, M. Kanoria / Eur. J. Pure Appl. Math, 4 (2011), 304-321 308 thermo-elasticity without energy dissipation (TEWOED(GN-II)) theory. The stress equation of motion in plain polar co-ordinates is given by ∂ τr r ∂ r + 1 r (τr r −τθθ ) = ρ ∂ 2u ∂ t2 . (4) For transversely isotropic body β11 = β22 [9]. Introducing the following dimensionless quan- tities U = C11 aβ11T0 u, (R,S) = � r a , b a � , (σR,σθ ) = 1 β11T0 (τr r ,τθθ ), Θ = T T0 , η = Gt a , (C1, C2) = 1 C11 (C11, C12), G2 = C11 ρ , (τ′q,τ′ν ,τ ′ T ) = G a � τq,τν ,τT � equations (1), (2), (3) and (4) become σR = ∂ U ∂ R + C2 U R −Θ, (5) σθ = C2 ∂ U ∂ R + U R −Θ, (6) � C2 T + (C 2 K +τ ′ νC 2 T ) ∂ ∂ η +τ′T C2 K ∂ 2 ∂ η2 �� ∂ 2Θ ∂ R2 + 1 R ∂Θ ∂ R � = � 1+τ′q ∂ ∂ η + 1 2 τ′2q ∂ 2 ∂ η2 �� ∂ 2Θ ∂ η2 + ε ∂ 2 ∂ η2 � ∂ U ∂ R + U R �� (7) and ∂ 2U ∂ R2 + 1 R ∂ U ∂ R − U R2 = ∂Θ ∂ R + ∂ 2U ∂ η2 , (8) where C2 T = K∗ ρCeG2 , C2 K = K aρCeG , and ε = β2 11T0 ρCeC11 are dimensionless constants, ε is the thermo- elastic coupling constant, CT is the nondimensional thermal wave velocity and CK is the damp- ing co-efficient. When CK = 0 and τT = τq = τν = 0, the corresponding equations become the equations of thermo-elasticity without energy dissipation theory (TEWOED(GN-II)) for Green-Naghdi model. The thermal and mechanical boundary conditions are given by qin = − ∂Θ ∂ R , U = 0 on R= 1, (9) Θ = 0, σR = 0 on R= S. (10) The term qin = qin(t) in equation (9) is the time dependent heat flux applied to the inner boundary (R= 1) of the disc. We assume that the medium is at rest and undisturbed initially. The initial conditions can be written as U = ∂ U ∂ η = ∂ 2U ∂ η2 = 0 and Θ = ∂Θ ∂ η = ∂ 2Θ ∂ η2 = 0 at η= 0,1≤ R≤ S. (11) A. Kar, M. Kanoria / Eur. J. Pure Appl. Math, 4 (2011), 304-321 309 3. Method of Solution Let ¦ U(R, p),Θ(R, p) © = ∫ ∞ 0 � U(R,η),Θ(R,η) e−pηdη (12) with Re(p) > 0 denote the Laplace transform of U and Θ respectively. On taking the Laplace transform, equations (7) and (8) reduce to d2Θ dR2 + 1 R dΘ dR = C � Θ+ ε � dU dR + U R �� (13) and d2U dR2 + 1 R dU dR − U R2 = p2U + dΘ dR , (14) where C = p2 � 1+τ′qp+ 1 2 τ ′2 q p2 � (1+τ′ν p)C2 T + p(1+τ′T p)C2 K . (15) 3.1. Eigen-value Approach Differentiating equation (13) with respect to R and using equation (14) we get d2 dR2 � dΘ dR � + 1 R d dR � dΘ dR � − 1 R2 � dΘ dR � = C � εp2U + (1+ ε) dΘ dR � . (16) From equations (14) and (16) we have the vector-matrix differential equation as follows LeV = eAeV , (17) where L ≡ d2 dR2 + 1 R d dR − 1 R2 , (18) and eV = � U , dΘ dR �T , eA= � D11 D12 D21 D22 � , (19) D11 = p2, D12 = 1, D21 = Cεp2 and D22 = C(1+ ε). To solve the vector matrix differential equation (17) we assume that eV = eX (m)ω(R, m), (20) where m is a scalar, eX is a vector independent of R and ω(R, m) is a non-trivial solution of the scalar differential equation Lω = m2ω. (21) A. Kar, M. Kanoria / Eur. J. Pure Appl. Math, 4 (2011), 304-321 310 The solution of the equation (21) is ω = [A1 I1(mR) +A2K1(mR)]. (22) Substituting equations (20) and (21) into the equation (17) we get eAeX = m2 eX , (23) where eX (m) is the eigenvector corresponding to the eigenvalue m2. The characteristic equation corresponding to eA can be written as m4 − (D11 + D22)m 2 + (D11D22 − D12D21) = 0. (24) The roots of the characteristic equation (24) are of the form m2 = m2 1 and m2 = m2 2, where m2 1 +m2 2 = D11 + D22, m2 1m2 2 = D11D22 − D12D21. (25) The eigenvectors X (m j), j = 1,2 corresponding to eigenvalues m2 j , j = 1,2 can be calculated as eX (m j) = � X1(m j) X2(m j) � = � D12 −(D11 −m2 j ) � , j = 1,2. (26) Therefore, the equation (20) can be written as eV = eX (m j)[A1 I1(m1R)+ B1K1(m1R)] + eX (m j)[A1I1(m2R)+ B1K1(m2R)], j = 1,2. (27) Therefore, from the equations in (19) we can write U = ∑ i=1,2 [Ai I1(miR) + BiK1(miR)] (28) and dΘ dR = − ∑ i=1,2 (p2 −m2 i )[Ai I1(miR) + BiK1(miR)], (29) where I1(miR) and K1(miR) are the modified Bessel functions of order one of first and second kind respectively. Ai and Bi ’s i = 1,2 are independent of R but depend on p and are to be determined from the boundary conditions. Using the recurrence relations of modified Bessel functions [34] we obtain from the equation (29) Θ = ∑ i=1,2 (p2 −m2 i ) mi [Ai I0(miR)+ BiK0(miR)], (30) since P1(mR) = − d dR � P0(mR) m � ; P = I , K . (31) A. Kar, M. Kanoria / Eur. J. Pure Appl. Math, 4 (2011), 304-321 311 Taking the Laplace transform of the equations (5) and (6) and using equations (28) and (30) we get σR = − ∑ i=1,2 Ai � 1− C2 R I1(miR)+ p2 mi I0(miR) � − ∑ i=1,2 Bi � 1− C2 R K1(miR) + p2 mi K0(miR) � (32) and σθ = ∑ i=1,2 Ai � 1− C2 R I1(miR) + ¨ (1− C2)mi − p2 mi « I0(miR) � + ∑ i=1,2 Bi � 1− C2 R K1(miR)+ ¨ (1− C2)mi − p2 mi « K0(miR) � . (33) Using the boundary conditions U = 0, dΘ dR = −qin on R = 1 and Θ = 0, σR = 0 on R = S and using the recurrence relations [34] from the equations (28), (29), (30) and (32) we obtain   A1 A2 B1 B2   =   W11 W12 W13 W14 W21 W22 W23 W24 W31 W32 W33 W34 W41 W42 W43 W44   −1  0 −qin 0 0   , (34) where W1i = P1(m j), W2i = −(p2−m2 j )P1(m j), (35) W3i = p2 −m2 j m j P0(m jS), W4i = 1− C2 S P1(m jS) + p2 m j P0(m jS) (36) where P = I for i = j =1,2; P = K for i = 3, j = 1 and i = 4, j = 2. From the equation (24) we obtain m1, m2 = 1 2 ( p α±pβ), where α,β = (p±pC)2+ Cε. (37) Therefore, m1 and m2 are real and positive quantities. 3.2. Finite element analysis In the finite element method, total domain is divided into a finite set of subintervals, i.e. line elements along the radial direction of the disc. The Galerkin finite element method is used to derive the stiffness and force matrices for the base element (e). For any base element (e) the displacement U and the temperature Θ can be approximated as U = [N](e){U∗}(e), Θ = [N](e){Θ∗}(e), (38) A. Kar, M. Kanoria / Eur. J. Pure Appl. Math, 4 (2011), 304-321 312 where [N](e) is the shape function approximating the displacement and temperature fields in the base element (e). The matrices {U∗}(e) and {Θ∗}(e) represent the nodal values of the displacement and temperature respectively. Using the equation (38) and the Galerkin finite element method over the volume of the base element V (e), equations (14) and (13) become ∫ V (e) �¨ d2 dR2 + 1 R d dR − 1 R2 − p2 « U − dΘ dR � NmdV = 0 (39) and ∫ V (e) �¨ d2 dR2 + 1 R d dR − C « Θ− εC � d dR + 1 R � U � NmdV = 0, (40) where Nm is the shape function. Equations (39) and (40) may be reduced to the weak form. Using the local coordinates R∗ = R− Ri, where Ri is the radius of the ith node of the base element, Equations (39) and (40) reduce to (dropping the asterisk for convenience) ∫ L 0 ¨ − �� 1 R+Ri d dR − 1 (R+Ri) 2 − p2 � U − dΘ dR � Nm(R+Ri) + d dR {(R+ Ri)Nm} dU dR « dR= (R+Ri)Nm dU dR ��� L 0 , (41) ∫ L 0 � − �� 1 R+ Ri d dR − C � Θ− εC � d dR + 1 R+ Ri � U � Nm(R+ Ri) + d dR {(R+ Ri)Nm} dΘ dR « dR= (R+Ri)Nm dΘ dR ��� L 0 , (42) where L = R j − Riis the length of the elements along the radial direction. The right-hand side terms of the equations (41) and (42) cancel each other between two adjacent elements, except the first node of the first element and the last node of the last element of the solution domain. These two nodes are located on the inner and outer boundaries of the disc. Thus applied boundary conditions may be considered as qin = − dΘ dR ��� 1 , U1 = 0, (43) ΘM = 0, S dU dR ��� M = −C2U M , (44) where U1 and U M are the radial displacements at the inner and outer boundaries of the disc respectively. The subscripts 1 and M denote the first and last nodes of the solution domain respectively. Substituting the equation (38) into the equations (41) and (42) we obtain (dropping the asterisks for convenience) ∫ L 0 Un � d dR {(R+ Ri)Nm} dNn dR − Nm dNn dR + NmNn R+Ri + p2(R+Ri)NmNn � dR A. Kar, M. Kanoria / Eur. J. Pure Appl. Math, 4 (2011), 304-321 313 + ∫ L 0 Θn � (R+Ri)Nm dNn dR � dR= (R+ Ri)Nm dU dR ��� L 0 , (45) ∫ L 0 Un � εC � (R+Ri)Nm dNn dR +NmNn �� dR+ ∫ L 0 Θn � d dR {(R+ Ri)Nm} dNn dR −Nm dNn dR + C(R+ Ri)NmNn � dR= (R+ Ri)Nm dΘ dR ��� L 0 , (46) Now the transfinite element equations (45) and (46) are expressed in the matrix form as � (K11) (K12) (K21) (K22) �� U Θ � = � F Q � . (47) The submatrices (K11), (K12), (K21), (K22), F , and Q are Kmn 11 = ∫ L 0 � d dR {(R+Ri)Nm} dNn dR − Nm dNn dR + NmNn R+ Ri + p2(R+ Ri)NmNn � dR, (48) Kmn 12 = ∫ L 0 � (R+ Ri)Nm dNn dR � dR, (49) Kmn 21 = ∫ L 0 εC � (R+Ri)Nm dNn dR + NmNn � dR, (50) Kmn 22 = ∫ L 0 � d dR {(R+ Ri)Nm} dNn dR − Nm dNn dR + C(R+Ri)NmNn � dR, (51) F =   0 0 · · 0 −C2UM   and Q =   qin 0 · · · 0   . (52) 4. Numerical Results and Discussions The solution in the space-time domain is obtained numerically by using Bellman et al. [4] method for fixed value of the space variable and for η = ηi , i = 1(1)7, where ηi ’s are com- puted from roots of the shifted Legendre polynomial of degree 7 (see Appendix). The compu- tations for the state variables are carried out for different values of R (1 ≤ R ≤ S) and values of ηi = 0.0257750,0.138382,0.352509,0.693147,1.21376,2.04612,3.67119. The material chosen for numerical evaluation is Sapphire. The physical data for Sapphire are† ρ = 3.96× 103kg/m3, ε= .0168, C11 = 4.96× 1010N/cm2, C12 = 1.15× 1010N/cm2 †as recorded on www.matweb. om A. Kar, M. Kanoria / Eur. J. Pure Appl. Math, 4 (2011), 304-321 314 and the relaxation time parameters and material parameters are taken as τq = 2.0×10−7sec, τT = 1.5× 10−7sec, τν = 1.0× 10−7sec, whereas CT = 2 and CK = 1.2. The dimensionless inner and outer radii are taken as R = 1 and S = 2 respectively. The di- mensionless input heat flux is defined as the Heaviside unit step function. That is qin(η) = 0 for η ≤ 0 and qin(η) = 1/s for η > 0. The results of the numerical evaluation of the thermo-elastic stresses, displacement and tem- perature are illustrated in figs. 1 to 8 using the Galerkin finite element method and eigen value approach. For the Galerkin finite element method the graphs are represented by FIN3P and FINGN-II respectively whereas for eigen-value approach these are denoted by 3PHASE and GN-II respectively. In all the figures 1 to 8 we observe that both the methods (eigen-value approach and Galerkin finite element method) produce similar qualitative behavior but the magnitudes differ slightly. This is obvious since the two methods are independent. Finite element method is fully numerical method, whereas, the eigen-value approach is partially analytical and partially numerical. Fig.1 shows the radial stress σR along the radial distance R for time η = .35. Since, by the assumption, the inner and outer boundaries of the disc are constrained and traction free respectively, the disc expand outwards. Also heat flux is applied on the inner boundary.Thus the stress waves propagate from inner boundary to outer boundary. It is also observed that the radial stress is compressive due to the application of the thermal shock on the inner boundary. Also the stress in each case is found to vanish on the outer boundary, which agrees with the imposed boundary condition. Fig.2 is plotted to show the variation of the circumferential stress σθ along R for fixed time η = .35. The elastic wave front is seen near the outer boundary.It is noticed that circumfer- ential stress in each case is compressive. Also each σθ attains its maximum magnitude on the inner boundary, since inner boundary is rigid and the heat flux is applied on it. Fig.3 is drawn to show the variation of the displacement U versus R for time η = .35. It is seen that the boundary condition for the displacement is satisfied on the inner boundary. The magnitudes of the displacements for the GN-II three-phase-lag models are almost the same in 1 ≤ R ≤ 1.2 whereas, the displacements propagate with different magnitudes in 1.2 < R ≤ 2. It is also noticeable that the magnitudes of the displacements for the case of the GN-II model asymptotically tend to zero much earlier than for those three-phase-lag model. Fig.4 depicts radial variation of the temperature Θ for fixed time η = .35. The boundary con- dition is satisfied on the outer boundary. It is also observed that the magnitude of each Θ is large in the case of three-phase-lag model in comparison to the GN-II model. Figs.5-6 show the stresses σR and σθ versus time η for fixed R = 1.4. It is noticed that when the time is small, the magnitudes of both the stresses for three-phase-lag model are large in comparison with those of the GN-II model. The opposite behavior is observed when time in- creases. This is due to the fact that there is no energy dissipation term in the GN-II model. It is also noticeable that as time increases both the stresses change their sign due to the reflection from the outer boundary. Fig.7 is plotted for displacement U against η for R= 1.4. The magnitude of each displacement is large for the GN-II model in comparison with those of three-phase-lag model. Fig.8 depicts the variation of temperature Θ versus time η for R= 1.4. Here the magnitude of the temperature for the GN-II model is less in 0 ≤ η ≤ .75 and 1.8 ≤ η ≤ 3.67 in comparison A. Kar, M. Kanoria / Eur. J. Pure Appl. Math, 4 (2011), 304-321 315 with those of three-phase-lag model, whereas, in the domain .75 < η < 1.8, the magnitude of the temperature for the GN-II model is large in comparison to that other model. For all the above numerical calculations FORTRAN − 77 programming Language has been used. Figure 1: Radial Stress vs R for η= 0.35. Figure 2: Cir umferential Stress vs R for η = 0.35. Figure 3: Radial Displa ement vs R for η = 0.35. A. Kar, M. Kanoria / Eur. J. Pure Appl. Math, 4 (2011), 304-321 316 Figure 4: Variation of Temperature vs R for η = 0.35. Figure 5: Radial Stress vs Time for R= 1.4. Figure 6: Cir umferential Stress vs Time for R= 1.4. REFERENCES 317 Figure 7: Displa ement vs Time for R= 1.4. Figure 8: Temperature vs Time for R= 1.4. ACKNOWLEDGEMENT. We are grateful to Prof. S.C. Bose of the Department of Applied Mathematics, University of Calcutta, for his kind help and guidance in the preparation of the paper. References [1] A Bagri, M Eslami. Generalized coupled themoelasticity of disks based on Lord-Shulman model, J. Thermal Stresses, 27 (2004) 691-704. [2] A Bagri, M Eslami. Analysis of thermoelastic waves in functionally graded hollow spheres based on the Green-Lindsay theory, J. Thermal Stresses, 30 (2007) 1175-1193. [3] A Bagri, M Eslami. A unifed generalized themoelasticity formulation; application to thick functionally graded cylinders, J. Thermal Stresses, 30, 911-930. 2007. [4] R Bellman, R Kolaba. J Lockette, Numerical inversion of the Laplace transform, American Elsevier Pub. Co., New York, 1966. REFERENCES 318 [5] D Chandrasekharaiah. Thermoelastic plane waves without energy dissipation, Mech. Res. Comm., 23 (1996) 549-555. [6] D Chandrasekharaiah. A note on the uniqueness of solution in the linear theory of ther- moelasticity without energy dissipation, J. Elasticity, 43 (1996) 279-283. [7] D Chandrasekharaiah. Hyperbolic thermoelasticity: A review of recent literature, Appl. Mech. Rev., 51 (1998), 8-16. [8] N Das, A Lahiri. Thermoelastic interactions due to prescribed pressure inside a spherical cavity in an unbounded medium, Int. J. Pure and Appl. Math. 31 (2000) 19-32. [9] R Dhaliwal, A Singh. Dynamic Coupled Thermoelasticity, Hindustan publication, 1980. [10] M Ghosh, M Kanoria. Generalized thermoelastic problem of a spherically isotropic elastic medium containing a spherical Cavity, J. Thermal Stresses, 31, 969-981. 2008. [11] M Ghosh, M Kanoria. Analysis of thermoelastic response in functionally graded spherically isotropic hollow sphere based on Green Lindsay theory, Acta Mech., 207, 51-67. 2009. [12] A Green, K Lindsay. Thermoelasticity, J. Elast. 2, 1-7.1972. [13] A Green, P Naghdi. A Re-examination of the basic postulate of thermo-mechanics, Proc. Roy. Soc. London, 432, 171-194. 1991. [14] A Green, P Naghdi. An unbounded heat wave in an elastic solid, J. Thermal Stresses, 15, 253-264. 1992. [15] A Green, P Naghdi. Thermoelasticity without energy dissipation, J. Elasticity, 31, 189-208. 1993. [16] R Hetnarski, J Ignaczak. Generalized thermoelasticity: Response of semi-space to a short laser pulse, J. Thermal Stresses, 17, 377-396. 1994. [17] A Kar, M Kanoria. Thermoelastic interaction with energy dissipation in an infinitely ex- tended thin plate containing a circular hole, Far East J. Appl. Math., 24, 201-217. 2006. [18] A Kar, M Kanoria. Thermoelastic interaction with energy dissipation in an unbounded body with a spherical hole, Int. J. Solids and Structures, 44, 2961-2971.2007. [19] A Kar, M Kanoria. Thermoelastic interaction with energy dissipation in a transversely isotropic thin circular disc, Euro. J. Mech. A/Solids, 26, 969-981.2007. [20] A Kar, M Kanoria. Generalized thermo-visco-elastic problem of a spherical shell with three- phase-lag effect, Appl. Math. Modelling, 33, 3287-3298. 2009. [21] A Lahiri, N Das and P Sen. Eigen value approach to three dimensional generalized ther- moelasticity, Bull. Cal. Math. Soc., 98, 305-318. 2005. REFERENCES 319 [22] S Lekhnitskii. Theory of elasticity of an anisotropic body, Mir publication, Moskow, 1980. [23] H Lord, Y Shulman. A generalized dynamic theory of thermoelasticity, J. Mech. Phys. Solids, 15, 299-309. 1967. [24] S Mallik, M Kanoria. Generalized thermo-elastic functionally graded sold with a periodi- cally varying heat source, Int. J. Solids and Stuctures, 44, 7633-7645. 2007. [25] Y Ootao, Y Tanigawa. Transient thermo-elastic problem of a functionally graded cylindrical pannel due to non-uniform heat supply, J. Thermal Stresses, 30, 441-457. 2007. [26] S Roychoudhuri. One-dimensional thermoelastic waves in elastic half-space with dual- phase-lag effects, J. of Mech. of Materials and Structures, 2, 489-503. 2007. [27] S Roychoudhuri. On a thermoelastic three-phase-lag model, J. Thermal Stresses, 30, 231- 238. 2007. [28] S Roychoudhuri and N Bandyopadhyay. Thermoelastic wave propagation in a rotating elastic medium without energy dissipation, Int. J. Math. and Math. Sci., 1, 99-107. 2004. [29] S Roychoudhuri, P Dutta. Thermoelastic interaction without energy dissipation in an infi- nite solid with distributed periodically varrying heat sources, Int. J. Solids and Structures, 42, 4192-4203. 2005. [30] Z Shao, T Wang, K Ang. Transient thermo-mechanical analysis of functionally graded hol- low cylinders, J. Thermal Stresses, 30. 81-104. 2007. [31] D Tzou. An unified field approach for heat conduction from macro to micro scales, ASME J. Heat Transfer, 117 8-16. 1995. [32] H Wang, H Ding, Y Chen. Thermo-elastic dynamic solution of a multilayered spherically isotropic hollow sphere for spherically symmetric problems, Acta Mech., 173, 131-145. 2005. [33] B Wang, Y Mai. Transient one dimensional heat conduction problems solved by finite ele- ment, Int. J. Mech. Sci., 47 303-317. 2005. [34] J Watson. Theory of Bessel functions(2nd edn.), Cambridge Univ. Press, 1980. Appendix Let the Laplace transform of σi(R,η) be given by σ j(R, p) = ∫ ∞ 0 e−pησ j(R,η)dη. (53) REFERENCES 320 We assume that σ j(R,η) is sufficiently smooth to permit the use of the approximate method we apply. Putting x = e−η in equation (53) we obtain σ j(R, p) = ∫ 1 0 x p−1 g j(R, x)d x , (54) where g j(R, x) = σ j(R,−log x). (55) Applying the Gaussian quadrature rule to the equation (54) we obtain the approximate rela- tion n∑ i=1 Wi x p−1 i g j(R, x i) = σ j(R, p), (56) where x i ’s(i = 1,2, . . . , n) are the roots of the shifted Legendre polynomial and Wi ’s (i = 1,2, . . . , n) are the corresponding weights [4] and p = 1(1)n. For p = 1(1)n, the equations (56) can be written as W1 g j(R, x1) +W2 g j(R, x2) + . . .+Wn g j(R, xn) = σ j(R, 1) W1 x1 g j(R, x1) +W2 x2 g j(R, x2) + . . .+Wn xn g j(R, xn) = σ j(R, 2) . . . . . . W1 xn−1 1 g j(R, x1) +W2 xn−1 2 g j(R, x2) + . . .+Wn xn−1 n g j(R, xn) = σ j(R, n) Therefore   g j(R, x1) g j(R, x2) · · · g j(R, xn)   =   W1 W2 · · · Wn W1 x1 W2 x2 · · · Wn xn · · · · · · · · · · · · · · · · · · W1 xn−1 1 W2 xn−1 2 · · · Wn xn−1 n   −1  σ j(R, 1) σ j(R, 2) · · · σ j(R, n)   . (57) (As the matrix is the product of diag {Wi}multiplied by Vander Monde matrix, it can be shown that the matrix is non-singular). Hence g j(R, x1), g j(R, x2), . . ., g j(R, xn) are known. From equations in (57) we can calculate the discrete values of g j(R, x i) i,e, σ j(R,ηi); (i =1,2,. . . ,7) and finally using interpolation we obtain the stress components σi(R,η); (i = R,θ). REFERENCES 321 Table 1: Roots and Weights of the Shifted Legendre Polynomial. n Roots Corresponding Weights 1 2.5446043828620886E-2 6.4742483084434816E-2 2 1.2923440720030282E-1 1.3985269574463828E-1 3 2.9707742431130145E-1 1.9091502525255938E-1 4 5.0000000000000000E-1 2.0897959183673466E-1 5 7.0292257568869853E-1 1.9091502525255938E-1 6 8.7076559279969706E-1 1.3985269574463828E-1 7 9.7455395617137909E-1 6.4742483084434816E-2