ARTICLE IN PRESS International Journal of Mechanical Sciences 49 (2007) 1325–1335 www.elsevier.com/locate/ijmecsci A unified generalized thermoelasticity; solution for cylinders and spheres A. Bagria, M.R. Eslamib, a Amirkabir University of Technology, P.O. Box 15875-4413 Tehran, Iran b Mechanical Engineering Department, Amirkabir University of Technology, Tehran, Iran Received 18 November 2006; received in revised form 29 March 2007; accepted 15 April 2007 Dedicated to Professor Franz Ziegler for his 70th birthday Available online 20 April 2007 Abstract A new unified formulation for the generalized theories of the coupled thermoelasticity based on the Lord–Shulman, Green–Lindsay, and Green–Naghdi models is proposed in this paper. The unified form of the governing equations is presented by introducing the unifier parameters. The formulations are derived and given for the anisotropic heterogeneous materials. The unified equations are reduced for the isotropic and homogeneous materials. Transforming the governing equations into the Laplace domain, they are analytically solved in the space domain for a hollow sphere and cylinder, where a parameter is introduced to consolidate the solution for the sphere and cylinder in a unified form. A thermal shock load is applied to the inner surface of the sphere and cylinder and the results are presented using a numerical inversion technique of the Laplace transform. The results are validated with the known data in the literature. r 2007 Elsevier Ltd. All rights reserved. Keywords: Generalized thermoelasticity; Unified formulation; Thick cylinders; Thick spheres; Thermal shock 1. Introduction is reported by Green and Lindsay [4]. In the Green–- Lindsay (GL) theory the Duhamel–Neumann relationships The conventional theory of thermoelasticity is based on and the entropy relation are modified by introducing two the Fourier’s heat conduction law. Due to the parabolic relaxation times that relate the stress and entropy to the nature of the energy equation of this theory, infinite temperature rate. An alternative approach in the formula- propagation speeds for the thermal disturbances are tion of a theory predicting the finite propagation speed of predicted. The concept of the hyperbolic nature involving the thermal disturbances is due to Green and Naghdi finite speeds of thermal disturbance is reported by Maxwell (GN), where they formulated three models of thermo- [1] for the first time, known as the second sound. Chester elasticity for homogeneous and isotropic materials [5,6] [2] provides some justification to the fact that the so-called labeled as models I, II, and III. These theories of second sound must exist in any solid. Most of the thermoelasticity (LS, GL, and GN theories) are known as approaches that came out to overcome the unacceptable the generalized theories, or thermoelasticity theories with prediction of the classical theory are based on the general the second sound effect or with finite thermal wave speed. notion of relaxing the heat flux in the classical Fourier heat Ignaczak [7] suggested a combined system of coupled conduction equation, thereby introducing a non-Fourier equations for the LS and GL theories. Also, the same effect. One of the simplest forms of these equation is due to author reported a survey of the domain of influence for the the work of Lord and Shulman [3]. In the Lord–Shulman results of the LS and GL theories [8]. Francis [9], Ignaczak (LS) theory a relaxation time is introduced and the [7], and Chandrasekharaiah [10,11] have reported brief Fourier’s heat conduction equation is modified. Another reviews of these theories. thermoelasticity theory that admits the second sound effect In situations such as those involving very short transient duration and/or when the sudden high heat flux is applied Corresponding author. Tel.: +98 21 6419736. to a structure, the second sound effect is important. Nayfeh E-mail address:
[email protected](M.R. Eslami). [12] and Nayfeh and Nemat-Nasser [13] used the Lord and 0020-7403/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmecsci.2007.04.004 ARTICLE IN PRESS 1326 A. Bagri, M.R. Eslami / International Journal of Mechanical Sciences 49 (2007) 1325–1335 Shulman theory to study the effects of the thermal coupling are transformed into the Laplace domain. Analytical on the plane harmonic thermoelastic waves in unbounded solution of the equations are obtained for hollow sphere media and Rayliegh surface waves propagating along the and cylinder, where a term is introduced to consolidate the free surface of a half space. Also, the important role of the solutions for the sphere and cylinder in a unified form. A coupling term in the response of elastic solids under thermal shock load is applied to the inner surface of the transient thermal loads is investigated by Amin and sphere and cylinder, where the solution in analytical form Sierakowski [14]. Chen and Dargush [15] used the is obtained in the space domain. The numerical inversion boundary element method to analyze the transient and of the Laplace transform is then employed to obtain the dynamic problems in generalized thermoelasticity of a half- solution in the time domain. space using the Laplace transform method. Chen and Lin [16] employed a hybrid numerical method based on the 2. Unified formulations of the LS, GL, and GN theories Laplace transform and control volume method to analyze the transient coupled thermoelastic problems with relaxa- The fundamental equations of the LS, GL, and GN tion times involving a nonlinear radiation boundary theories in a unified form are presented here introducing condition. Boundary element method is employed by the terms Z and t3 as unifier parameters. These equations in Hosseini Tehrani and Eslami [17,18] for the analysis of general form are: coupled thermoelastic problems in a finite domain, where Equations of motion: they studied the coupling coefficient effects on the thermal r r þ rb ¼ r€u, (1) and elastic waves propagation. A transfinite element method is considered by Bagri and Eslami [19] to study Linear strain–displacement relations: the generalized coupled thermoelasticity of an annular disk E ¼ 12ðru þ ðruÞ0 Þ, (2) based on the LS theory. They investigated the thermal and stress waves propagation through the radius of the disk Hooke’s law for the linear thermoelastic materials: and showed that for thermal shock problems the coupling _ r ¼ C E bðT T 0 þ t1 TÞ, (3) coefficient has significant effect on the variation of thermal stresses, displacement, and temperature. Energy balance equation: Several investigators employed the GN models to solve a _ r q ¼ R T 0 S, (4) variety of thermoelastic problems. The uniqueness of solution of the governing equations for the GN theory Entropy relationship: formulated in terms of stress and energy-flux is established rc 1 b in Ref. [20]. Chandrasekharaiah [21] studied the one- S¼ ðT þ t2 T_ T 0 Þ þ b : E C rT, (5) T0 T0 dimensional thermal wave propagation in a half-space based on the GN model due to a sudden exposure of Heat conduction equation: temperature to the boundary, using the Laplace transform b T, Z q þ Z s q_ þ t3 q_ ¼ Z K rT t3 K rT_ t3 K rT C _ method. Chandrasekharaiah [22] also presented the com- plete solutions of the governing field equations for the GN (6) theory. Sharma and Chauhan [23] investigated the dis- where r is the mass density, r is the Cauchy’s stress tensor, turbances produced in a half-space under the application of u is the displacement vector, b is the body force vector per a point load and thermal source acting on the boundary of unit mass, q is the heat flux vector, T 0 is the reference the half-space. The material is assumed to be homogeneous temperature, T is the absolute temperature, S is entropy and isotropic. The Laplace and Hankel transforms are used per unit volume, R is the internal heat source per unit and different theories of generalized thermoelasticity are volume per unit time, E is the strain tensor, b is the second employed to provide a basis to compare the results. Taheri order tensor of stress–temperature moduli, K is the second et al. [24] studied the problem of coupled thermoelasticity order tensor of thermal conductivity, C is the forth order of a layer based on the GN theory. The problem was tensor of elastic moduli, and c is the specific heat, transformed into the Laplace domain, where the analytical respectively. Also, s is the second order tensor of relaxation solution was obtained. An inverse numerical method was times in the LS model, t1 and t2 are the relaxation times then employed to obtain the solution in real time domain. and C b is a vector of new material constants proposed by In this paper, a new unified formulation for the Green and Lindsay, and K is the second order tensor of generalized coupled thermoelasticity theories based on new material constants associated with the GN theory. the LS, GL, and GN models is proposed. The unifier Also, Z and t3 are terms introduced to consolidate all parameters are introducing to consolidate the equations of theories into a unified system of equations. In Eqs. (1)–(6) the LS, GL, and GN theories into a single system of the superscript dot (.) denotes differentiation with respect equations for the anisotropic and heterogeneous materials. to time. Meanwhile, r is the del operator and indicates the The equations are also simplified for the isotropic and gradient of a function, ðr:Þ denotes the divergence homogeneous materials. The reduced equations for the operator, and the superscript prime ð0 Þ indicates the isotropic and homogeneous material are considered and transpose of the matrix. In Eqs. (1)–(6) the double dot ARTICLE IN PRESS A. Bagri, M.R. Eslami / International Journal of Mechanical Sciences 49 (2007) 1325–1335 1327 product of two second-order tensors means the trace of or in index notation their product. For example, b : E in index form is bij E ij . C ijkl ¼ ldij dkl þ mðdik djl þ dil djk Þ The product of CE is a second-order tensor and in index notation is C ijkl E kl , the product ofsq is a first-order tensor bij ¼ bdij ; tij ¼ t0 dij , and in index notation is tij qj , and C b rT is a scaler which K ij ¼ kdij ; K ij ¼ kdij , b in index notation is C i T ;i . Bij ¼ expðt=t0 Þdij ; eij ¼ expðe B t=t0 Þdij , Eqs. (1)–(6) include the governing coupled system of equations for the classical theory of thermoelasticity when bi ¼ 0, C s ¼ 0, Cb ¼ 0, t1 ¼ t2 ¼ t3 ¼ 0 and Z ¼ 1. The given system of equations yields the equations of the LS theory when where l and m are the Lame´ constants, k is the thermal Z ¼ 1, Cb ¼ 0 and t1 ¼ t2 ¼ t3 ¼ 0. The governing equations conductivity, t0 is the relaxation time proposed by Lord for the GL theory are obtained when s ¼ 0, Z ¼ 1 and and Shulman, indicates the tensor product, and dij is the t3 ¼ 0. To obtain the governing equations for the GN theory Kronecker delta. Eqs. (7) and (8) for the isotropic type III, we set s ¼ 0, C b ¼ 0, Z ¼ t1 ¼ t2 ¼ 0, and t3 ¼ 1. homogeneous medium appear in the form Two special cases of the GN theory, namely type II and I, Z t t3 kr:ðrTÞ_ þ t3 k r:ðrTÞ þ Z k expððt e tÞ=t0 Þr:ðrTÞ de t may be obtained from the equations of the GN theory type t0 0 III by setting K ! 0 and K ! 0, respectively. To obtain the rcðt2 þ t3 ÞT€ ZrcT_ ZT 0 bðr:_uÞ GN theory type II, from the equations of the GN theory type III, we set K ! 0. When K ! 0 the equations of the GN t3 T 0 bðr:€uÞ þ ZR þ t3 R_ ¼ 0, ð11Þ theory type III reduce to the GN theory type I, which is 0 _ lrðr:uÞ þ mr:ððruÞ þ ðruÞ Þ brðT þ t1 TÞ þ rb ¼ r€u. ð12Þ identical with the classical theory of thermoelasticity. Also, these equations may be written as To obtain the system of equations in terms of displace- Z ment and temperature, terms r; E; q, and S may be 2 2 _ k t eliminated with the help of Eqs. (1)–(6). Considering the t3 k r T þ t3 kr T þ Z expððt e tÞ=t0 Þðr2 TÞ de t t0 0 values for the terms Z; t3 for each theory and noting that rcðt2 þ t3 ÞT€ ZrcT_ ZT 0 bðr:_uÞ t3 T 0 bðr:€uÞ the term s is only appeared in the LS model, the system of equations in terms of the displacement components and þ ZR þ t3 R_ ¼ 0, ð13Þ temperature are (see also Appendix A for more details) 2 _ þ rb ¼ r€u, mr u þ ðl þ mÞrðr:uÞ brðT þ t1 TÞ ð14Þ Z t _ t3 r:ðKrTÞ þ t3 r:ðK rTÞ þ Zr B 1 e Bs KðrTÞ de 1 t where r2 is the Laplacian operator. Moreover, the heat 0 conduction equation and the stress–displacement relations r c ðt2 þ t3 ÞT€ Z r c T_ Z T 0 b : ðru_ Þ t3 T 0 b : ðru€ Þ for the isotropic and homogeneous materials are þ 2C b ðrTÞ_ þ ðr CÞb T_ þ ZR þ t3 R_ ¼ 0, ð7Þ _ r ¼ lðr:uÞI þ mððruÞ þ ðruÞ0 Þ bðT T 0 þ t1 TÞI, ð15Þ r ðC ðr uÞÞ ðr bÞ T T 0 þ t1 T _ Zq þ Zt0 q_ þ t3 q_ ¼ ZkrT t3 krT_ t3 k rT. ð16Þ b r T þ t1 T_ þ rb ¼ r€u, ð8Þ where X1 n 3. Governing equations for symmetrically loaded spheres t 1 n B ¼ expðt s1 Þ ¼ ðs Þ , and cylinders n¼0 n! 1 en X Consider a symmetrically loaded sphere, or a long e ¼ expðe t B t s1 Þ ¼ ðs1 Þn . ð9Þ cylinder under plane strain condition, exposed to a n¼0 n! symmetric thermal shock load. Eqs. (13)–(16) may be Here, ðs1 Þ0 ¼ I, and I is the Identity tensor. written in spherical or cylindrical polar coordinates and in Reconsideration of Eqs. (7) and (8) for the GN theory dimensionless form. In order to obtain the equations in type II reveals that no damping term is appeared in the dimensionless form, the following dimensionless para- system of equations and therefore the GN theory type II is meters are introduced: known as the thermoelasticity without energy dissipation. For the isotropic materials, the material constants are 1 v 1 ðl þ 2mÞ x¯ ¼ x; ¯t ¼ t; u¯ ¼ u, ‘ ‘ ‘ bT d C ¼ lI I þ mðII þ I IÞ, v v v v ¯t0 ¼ t0 ; ¯t1 ¼ t1 ; ¯t2 ¼ t2 ; ¯t3 ¼ t3 , b ¼ bI; s ¼ t0 I, ‘ ‘ ‘ ‘ K ¼ kI; K ¼ k I, T T 0 1 1 T¯ ¼ ; r¯ ¼ r; q¯ ¼ q, Td bT d rcvT d e ¼ expðe B ¼ expðt=t0 ÞI; B t=t0 ÞI, b ¼ 0, R ¯ ¼ ‘ R; b¯ ¼ b‘T d b, ð17Þ C ð10Þ rcvT d ðl þ 2mÞv2 ARTICLE IN PRESS 1328 A. Bagri, M.R. Eslami / International Journal of Mechanical Sciences 49 (2007) 1325–1335 where ‘ is a characteristic length, v is a characteristic speed where u; br ; qr , and r are, respectively, the components of and T d is a characteristic temperature. Eqs. (13)–(16) in displacement, body force, heat flux, and coordinate terms of the dimensionless parameters, dropping the bar variable in the radial direction and srr and syy are the for convenience, take the form radial and hoop stresses. Also, term m is a parameter Z introduced to consolidate the equations in the polar and c2 t spherical coordinates in a unified form. When m ¼ 1 the t3 c2T r2 T þ t3 c2K r2 T_ þ Z K expððt etÞ=t0 Þðr2 TÞ de t t0 0 equations in polar coordinates, and when m ¼ 2 the ðt2 þ t3 ÞT€ ZT_ ½t3 r:€u þ Zr:_u þ ZR þ t3 R_ ¼ 0,ð18Þ equations in spherical coordinates are obtained, respec- _ þ b ¼ u€ , ð19Þ tively. c22 r2 u þ ðc21 c22 Þrðr:uÞ c21 ðrT þ t1 rTÞ To solve the system of Eqs. (23)–(26) the Laplace transform method may be used. Suppose that the body is initially at rest and at reference temperature and the initial c22 c2 _ r¼ 1 2 2 ðr:uÞI þ 22 ðru þ ðruÞ0 Þ ðT þ t1 TÞI, ð20Þ value of the heat flux and rate of temperature are zero. c1 c1 Thus uðr; 0Þ ¼ 0; _ 0Þ ¼ 0, uðr; Zq þ Zt0 q_ þ t3 q_ ¼ Zc2K rT t3 c2K rT_ t3 c2T rT, (21) _ 0Þ ¼ 0; Tðr; 0Þ ¼ 0; Tðr; qðr; 0Þ ¼ 0. ð27Þ where Applying the Laplace transform to Eqs. (23)–(26) leads ðl þ 2mÞ m k to c21 ¼ ; c22 ¼ 2; c2T ¼ , rv2 rv rcv2 2 2 k b2 T 0 q mq Zc2K q mq c2K ¼ ; ¼ . ð22Þ ðt3 c2T þ st3 c2K Þ 2 þ T þ þ T rc‘v rcðl þ 2mÞ qr r qr t0 s þ 1 qr2 r qr q 1 In Eqs. (18)–(21), c1 and c2 represent the nondimensional ðt2 þ t3 Þs2 T ZsT ðt3 s2 þ ZsÞ þ m u qr r speeds of purely elastic dilatational and shear waves, þ ðZ þ t3 sÞR ¼ 0, ð28Þ respectively. Moreover, cT and cK are, respectively, the speed of thermal wave disturbancespand the damping ffiffiffiffiffiffiffiffiffiffiffiffiffi coefficient in the GN model, while cK = t0 þ t2 is the speed q q 1 q of thermal wave disturbances in the LS and GL models and c21 2 þ m u c1 ð1 þ t1 sÞ T þ br ¼ s2 u , (29) qr qr r qr is the usual thermoelastic coupling parameter. Eqs. (18)–(21) in the spherical and cylindrical polar coordinates are 2 3 " # qu srr c22 c22 2 6 qr 7 q2 m q mq _ ¼2 4 5þ 12 2 2 q syy c21 c1 t3 c2T þ T þ t c 3 K þ T u =r qr2 r qr qr2 r qr Z 2 " 1 # c2K t q mq qu u þZ e expððt t Þ=t0 Þ þ T de t þm ð1 þ t1 sÞT , ð30Þ t0 0 qr2 r qr qr r 1 € _ q 1 q 1 ðt2 þ t3 ÞT ZT t3 þ m u€ þ Z þ m u_ qr r qr r qT _ þ ZR þ t3 R ¼ 0, ð23Þ ðZ þ Zt0 s þ t3 sÞqr ¼ ðZc2K þ t3 c2T þ st3 c2K Þ , (31) qr q q 1 qT qT_ c21 þ m u c21 þ t1 € þ br ¼ u, ð24Þ where s is the Laplace parameter and the asterisk sign qr qr r qr qr indicates the Laplace transform of the terms. Considering 23 the prescribed values for t0 ; t2 ; t3 ; and Z defined to obtain " # qu the equations for the LS, GL, and GN models, in the srr c22 6 7 c2 ¼ 2 4 qr 5 þ 1 2 22 absence of body force and heat source, Eqs. (28)–(31) may syy c21 c1 also be written as u=r " 1 # qu u 1 q m qT 2 1 qðrm u Þ þm _ ðT þ t1 TÞ , ð25Þ r o2 T ¼ z m , (32) qr r 1 rm qr qr r qr qT qT_ q 1 qðrm u Þ 2 qT Zqr þ Zt0 q_ r þ t3 q_ r ¼ ðZc2K þ t3 c2T Þ t3 c2K , (26) o1 u ¼ g , (33) qr qr qr rm qr qr ARTICLE IN PRESS A. Bagri, M.R. Eslami / International Journal of Mechanical Sciences 49 (2007) 1325–1335 1329 23 " # qu u ¼ ue =rðm1Þ=2 . These solutions are srr c22 c22 6 qr 7 X 2 ¼2 4 5þ 12 2 syy c21 c1 u ¼ rðm1Þ=2 fAi I ðmþ1Þ=2 ðDi rÞ þ Bi K ðmþ1Þ=2 ðDi rÞg, ð43Þ u =r i¼1 1# " qu u X 2 þm gT , ð34Þ T ¼ rðm1Þ=2 fG i I ðm1Þ=2 ðDi rÞ þ H i K ðm1Þ=2 ðDi rÞg, ð44Þ qr r 1 i¼1 where D2i ði ¼ 1; 2Þ are the roots of the characteristic 1 qT equation qr ¼ 2 , (35) o3 qr D4 ðo21 þ o22 þ x2 ÞD2 þ o21 o22 ¼ 0, (45) where and I ðmþ1Þ=2 and K ðmþ1Þ=2 are the modified Bessel functions of the first kind and second kind and order ðm þ 1Þ=2, s2 ðt0 þ t2 þ t3 Þs2 þ Zs respectively. Also, I ðm1Þ=2 and K ðm1Þ=2 are, respectively, o21 ¼ ; o22 ¼ , c21 ðZc2K þ t3 c2T þ st3 c2K Þ the modified Bessel functions of the first kind and second x2 ½ðt0 þ t3 Þs2 þ Zsg kind and order ðm 1Þ=2. In Eqs. (43) and (44) Ai ; Bi ; G i , z¼ ; x2 ¼ ; g ¼ 1 þ t1 s, and H i are the constants of integration that must be g ðZc2K þ t3 c2T þ st3 c2K Þ obtained using the given boundary conditions. ðZ þ Zt0 s þ t3 sÞ o23 ¼ . ð36Þ Now, consider the following relations for the derivatives ðZc2K þ t3 c2T þ st3 c2K Þ of the modified Bessel functions: Eliminating u in Eqs. (32) and (33) results in d½I n ðDi rÞ n ¼ Di I n1 ðDi rÞ I n ðDi rÞ dr r 1 q m q 1 q m q 2 n r r o2 T ¼ Di I nþ1 ðDi rÞ þ I n ðDi rÞ, rm qr qr rm qr qr r 2 1 q m q 2 d½K n ðDi rÞ n o1 m r o2 T ¼ Di K n1 ðDi rÞ K n ðDi rÞ r qr qr dr r n 2 1 q m qT ¼ Di K nþ1 ðDi rÞ þ K n ðDi rÞ. ð46Þ ¼x m r . ð37Þ r r qr qr Substitution of Eqs. (43) and (44) into Eq. (32), or (33), and Then considering Eqs. (46) results in ðr2 o21 Þðr2 o22 ÞT ¼ x2 r2 T , (38) Ai ðD2i o21 Þ Gi ¼ , ð47Þ gDi where Bi ðD2i o21 Þ Hi ¼ . ð48Þ 1 q m q q2 m q gDi 2 r ¼ m r ¼ 2þ . (39) r qr qr qr r qr Thus Also, eliminating term T in Eqs. (32) and (33) leads to the X 2 u ¼ rðm1Þ=2 fAi I ðmþ1Þ=2 ðDi rÞ þ Bi K ðmþ1Þ=2 ðDi rÞg, following equation in terms of u : i¼1 (49) q 1 q m q 1 qðrm u Þ 2 r o1 u qr rm qr qr rm qr 2 2 X D o2 2 q 1 qðrm u Þ 2 2 q 1 qðrm u Þ T ¼ rðm1Þ=2 i 1 ðAi I ðm1Þ=2 ðDi rÞ o2 o1 u ¼ x . gDi qr rm qr qr rm qr i¼1 ð40Þ Bi K ðm1Þ=2 ðDi rÞÞ . ð50Þ With definition of the operator n2 as Substituting Eqs. (49) and (50) into Eqs. (34) and (35), considering Eqs. (46) for the derivatives of the modified 2 q 1 q m q2 m q m n ¼ ðr Þ ¼ þ , (41) Bessel functions, the radial and hoop stresses and heat flux qr rm qr qr2 r qr r2 are obtained as Eq. (40) appears in the form X2 2 o1 srr ¼ rðm1Þ=2 ðAi I ðm1Þ=2 ðDi rÞ Bi K ðm1Þ=2 ðDi rÞÞ ðn2 o21 Þðn2 o22 Þu ¼ x2 n2 u . (42) i¼1 Di Solution of Eqs. (38) and (42) are obtained with c¯ 22 1 2m 2 ðAi I ðmþ1=2Þ ðDi rÞ þ Bi K ðmþ1Þ=2 ðDi rÞÞ , ð51Þ interchanging the variables as T ¼ Te =rðm1Þ=2 and c¯ 1 r ARTICLE IN PRESS 1330 A. Bagri, M.R. Eslami / International Journal of Mechanical Sciences 49 (2007) 1325–1335 X2 producing a thermal shock load applied to the structure. c¯ 2 1 syy ¼ rðm1Þ=2 2 22 ðAi I ðmþ1Þ=2 ðDi rÞ To carry out the numerical analysis, the following material i¼1 c¯ 1 r constants and parameters are considered: c¯ 2 o2 þ Bi K ðmþ1Þ=2 ðDi rÞÞ 2Di 22 1 c1 ¼ 1; c2 ¼ 0:535; ¼ 0:02. (55) c¯ 1 Di ðAi I m1 2 ðD i rÞ B K i ðm1Þ=2 ðD i rÞÞ , ð52Þ 4.1. Response of cylinder to heat flux shock load X 2 D2i o21 qr ¼ rðm1Þ=2 ðAi I ðmþ1Þ=2 ðDi rÞ Consider a thick long cylinder under a heat flux shock go23 i¼1 load qr ¼ f ðtÞ applied to its inner surface, where the þ Bi K ðmþ1Þ=2 ðDi rÞÞ. ð53Þ temperature change at the outer surface of the cylinder is assumed to be zero. Therefore, the boundary conditions in While the solution for the cylindrical polar coordinates the Laplace domain are and the spherical coordinates are presented in a unified form and in terms of the Bessel functions, solution for 10 000 qr ¼ ; u ¼ 0 at r ¼ a, the spherical coordinates may be given in terms of the sðs þ 100Þ2 exponential functions [25] and the solution for the cylindrical polar coordinates in terms of the Kelvin T ¼ 0; srr ¼ 0 at r ¼ b, (56) functions [26]. The Kelvin functions are introduced in where a and b are the inner and outer radii of the cylinder, terms of the cosine and sine functions. respectively. With the given boundary conditions and Now, the boundary conditions may be applied to obtain setting m ¼ 1, Eqs. (49)–(53) are used to obtain the the constants Ai and Bi . Eqs. (49)–(53) are in the Laplace constants Ai and Bi for the cylinder in the Laplace domain. To transform the equations into the physical time transform domain. The numerical inversion of the Laplace domain, the inverse Laplace transform should be used. The transform, proposed by Honig and Hirdes [27], may be inversion of the Laplace transform may be obtained employed to obtain the solutions in physical time domain. numerically, using the technique given by Honig and As an example, consider a cylinder with inner and outer Hirdes [27]. radii a ¼ 1 and b ¼ 2. For the LS model, the characteristic length l and characteristic speed v in Eq. (17) are chosen 4. Numerical results and discussions such that t0 ¼ 4 and cK ¼ 1. With the assumed numerical pffiffiffiffi values, c1 and cK = t0 , the speeds of propagation of elastic In this section, the dynamic responses of cylinder and and thermal disturbances, become unity and 0:5, respec- sphere are obtained using the derived formulations and tively. Figs. 2–4 show the distribution of temperature, solutions. To study the responses under different boundary radial, and hoop stresses along the radial direction at conditions, two types of thermal loads are considered. For different times. In the temperature distribution, times t ¼ the cylinder and sphere the mechanical boundary condi- 0:2 to t ¼ 1:8 show the propagation of the wave, while time tions are assumed to be fixed and traction free at the inner t ¼ 2:2 shows the reflection of the wave from the outer and outer surfaces, respectively. Also, the thermal input is boundary of the cylinder. In the stress distributions, times assumed to be a given function of nondimensional time as t ¼ 0:2 to t ¼ 1 indicate the elastic wave propagations and f ðtÞ ¼ 1 ð1 þ 100tÞe100t . (54) time t ¼ 1:3–1.8 depicts the reflection of elastic waves from the outer boundary. Figs. 2–4 obviously show the thermal The variation of function f ðtÞ versus time is shown in and elastic wave fronts. It may be seen from Fig. 2 that for Fig. 1. It is observed that f ðtÞ abruptly vary from zero to 1, Fig. 2. Distribution of nondimensional temperature along the radial Fig. 1. Input thermal load versus nondimensional time. direction of the cylinder for different times considering LS theory. ARTICLE IN PRESS A. Bagri, M.R. Eslami / International Journal of Mechanical Sciences 49 (2007) 1325–1335 1331 Fig. 3. Distribution of nondimensional radial stress along the radial Fig. 5. Distribution of nondimensional temperature along the radial direction of the cylinder for different times considering LS theory. direction of the cylinder for different times considering GL theory. Fig. 4. Distribution of nondimensional hoop stress along the radial Fig. 6. Distribution of nondimensional radial stress along the radial direction of the cylinder for different times considering LS theory. direction of the cylinder for different times considering GL theory. the LS theory temperature wave has a steep jump at thermal wave front. Moreover, Fig. 2 depicts that the thermal wave hits the position of r ¼ 1:1 at time t ¼ 0:2, which means that the thermal wave propagates with the speed of 0:5. Also, Figs. 3 and 4 show that the elastic wave hits the position of r ¼ 1:2 at time t ¼ 0:2, which means that the elastic wave propagates with the speed of unity. Meanwhile, it may be seen form Figs. 3 and 4 that the stress waves have two steep jumps at elastic and thermal waves fronts. As may be seen from Fig. 3, the sign of radial stress wave is reversed when the elastic wave reflection from the free boundary, i.e. the outer boundary, reaches Fig. 7. Distribution of nondimensional hoop stress along the radial the thermal wave front. This situation takes place at a time direction of the cylinder for different times considering GL theory. between t ¼ 1:3 and 1:4. For the GL theory, the characteristic length l and the characteristic speed v in Eq. (17) are chosen such that t1 ¼ GL theory and under the heat flux shock load no steep t2 ¼ 4 and cK ¼ 1. With the assumed numerical values for jump is revealed for the temperature wave. However, cK and the relaxation times (when c1 ¼ 1) the speed of similar stress wave propagations are obtained. Also, two propagation ofpffiffiffiffithe thermal disturbances for the GL theory steep jumps are observed for the stress waves. becomes cK = t2 ¼ 0:5. Figs. 5–7 show the distribution of The characteristic length l and characteristic speed v are temperature, radial, and hoop stresses along the radial chosen such that t3 ¼ 1 and cT ¼ 0:5 for the GN theory direction at different times for the GL theory. Interpreta- type II (cK ¼ 0). The GN theory type II reveals no tion of the speed of wave propagation and times of damping term and is suitable for analysis of the primary propagation and reflection of the waves in the GL theory stages of thermal shock applications in a structure. While are similar to those discussed for the LS theory. It may be the solutions are given for all types of the GN theory, the seen from Fig. 5 that in contrast to the LS theory, in the results are only shown for type II. Figs. 8–10 show the ARTICLE IN PRESS 1332 A. Bagri, M.R. Eslami / International Journal of Mechanical Sciences 49 (2007) 1325–1335 theory type II predicts the wave propagation without energy dissipation, thus the heat flux for larger values of radius is decreased causing decrease in the temperature and stress wave fronts. Reason for the stress and temperature amplitude reduction along the radial direction is that the circumference of the cylinder increases as the radius is increased resulting into lower values of heat flux. Comparison between the generalized coupled thermo- elasticity theories shows that under heat flux shock loads, the LS, GL, and GN theories predict similar patterns for the stresses, but different values for thermal filed. Also, it is seen from the figures that the LS and GN theories predict Fig. 8. Distribution of nondimensional temperature along the radial similar temperature values, which are larger than those direction of the cylinder for different times considering GN theory. obtained using the GL theory when the cylinder is exposed to a heat flux shock load. 4.2. Response of sphere to temperature shock load Consider a thick sphere under spherically symmetric thermal shock load. The behavior of the sphere based on the LS, GL, and the GN theories is investigated. A temperature shock load with the function T ¼ f ðtÞ is applied to the inner surface of the sphere, where the outer surface is insulated. Therefore, the boundary conditions in the Laplace domain are 10 000 T ¼ ; u ¼ 0 at r ¼ a, sðs þ 100Þ2 Fig. 9. Distribution of nondimensional radial stress along the radial direction of the cylinder for different times considering GN theory. qT ¼ 0; srr ¼ 0 at r ¼ b, (57) qr where a and b are the inner and outer radii of the sphere, respectively. The assumed boundary conditions are applied to Eqs. (49)–(53) to obtain the constant Ai and Bi . In this case, we set m ¼ 2 and employ the numerical inversion of the Laplace transform to obtain the temperature and stress fields in time domain. The numerical values for the inner and outer radii of the sphere are assumed as a ¼ 1 and b ¼ 2, respectively. Interpretations of the speed of wave propagations and reflections from the boundaries in the sphere based on the LS, GL, and GN theories, shown in Figs. 11–19, are similar Fig. 10. Distribution of nondimensional hoop stress along the radial direction of the cylinder for different times considering GN theory. thermal and elastic waves propagation along the radial direction. It may be observed from Fig. 8 that the wave front at t ¼ 0:2 hits the position of r ¼ 1:1 and the temperature wave has a steep jump at thermal wave front. Thus, the speed of the thermal wave propagation, i.e. cT ¼ 0:5, may be obtained from the results. Also, Figs. 9 and 10 reveal a speed of unity for the elastic wave propagation, where two abrupt variations are seen from Fig. 11. Distribution of nondimensional temperature along the radial the figures at thermal and elastic waves fronts. The GN direction of the sphere for different times considering LS theory. ARTICLE IN PRESS A. Bagri, M.R. Eslami / International Journal of Mechanical Sciences 49 (2007) 1325–1335 1333 Fig. 15. Distribution of nondimensional radial stress along the radial Fig. 12. Distribution of nondimensional radial stress along the radial direction of the sphere for different times considering GL theory. direction of the sphere for different times considering LS theory. Fig. 16. Distribution of nondimensional hoop stress along the radial Fig. 13. Distribution of nondimensional hoop stress along the radial direction of the sphere for different times considering GL theory. direction of the sphere for different times considering LS theory. Fig. 14. Distribution of nondimensional temperature along the radial Fig. 17. Distribution of nondimensional temperature along the radial direction of the cylinder for different times considering GL theory. direction of the sphere for different times considering GN theory. thermal and elastic wave fronts. The prior Dirac-delta-type to those given for the waves propagation in the cylinder. behavior of the solutions for the radial stress is in tension An interesting result obtained from Figs. 11–13 and Figs. and the latter is in compression, while both of them are 17–19 is that a similar pattern is obtained for the LS and compressive for the hoop stress. Moreover, it may be seen GN theories and therefore the LS and GN models reveal from the figures that for the LS and GN theories the radial close results in case where a temperature shock load is stress wave is changed from compression to tension when applied to the sphere. Considering the results obtained the reflection of the elastic wave front from the free from the GL theory, Figs. 14–16 indicate that when a boundary of the sphere reaches to the thermal wave front. temperature shock load is imposed to the sphere, the Moreover, for the GL theory Fig. 15 depicts that the sign temperature wave has a steep jump at the thermal wave of Dirac-delta-type behaviour of the radial stress wave at front. Also, the steep jump in the temperature wave front the elastic wave front is reversed from compression to results in a Dirac-delta-type solution for the stresses at tension when it is reached to the free boundary of the ARTICLE IN PRESS 1334 A. Bagri, M.R. Eslami / International Journal of Mechanical Sciences 49 (2007) 1325–1335 into the Laplace transform domain, they are analytically solved in the space domain for a hollow sphere and cylinder, where a parameter is introduced to consolidate the solutions for the sphere and cylinder in a unified form. Two types of thermal shock loads are applied to the inner surface of the sphere and cylinder, where the solution in analytical form is obtained in the space domain. The numerical inversion of the Laplace transform is then employed to obtain the solutions in the time domain. Results obtained from the figures show that under the heat flux loads the LS, GL, and GN theories predict similar patterns for the stresses with two steep jumps in the stress Fig. 18. Distribution of nondimensional radial stress along the radial distributions at thermal and elastic wave fronts. On the direction of the sphere for different times considering GN theory. other hand, they predict different values for the thermal filed. Also, it may be seen from the figures that the LS and GN theories predict similar temperature values, which are larger than those obtained using the GL theory. Compar- ison between the LS, GL, and GN theories shows that under the temperature shock loads, the LS and GN theories predict similar patterns for the temperature and stresses. Moreover, it may be found from the figures that the GL theory predicts a Dirac-delta-type result for the stresses due to the existence of term T_ in the stress–strain relations. Also, the temperature patterns obtained from these theories are similar. Appendix A. Derivation of Eq. (7) Fig. 19. Distribution of nondimensional hoop stress along the radial Applying del operator with dot sense, ðr : Þ, to Eq. (6) direction of the sphere for different times considering GN theory. results in Z r : q þ Z r : ðs q_ Þ þ t3 r : q_ sphere. On the other hand, in the GL theory the radial b TÞ. _ t3 r : ðK rTÞ r : ðC ¼ Z r : ðK rTÞ t3 r : ðK rTÞ _ stress wave sign is not reversed when the elastic wave reflection from the free boundary of the sphere has reached ð58Þ the thermal wave front. Substituting the first and third terms on the left-hand side Comparison between the LS, GL, and GN theories of the above equation by Eq. (4) gives shows that under the temperature shock loads, the LS and _ þ Z r : ðs q_ Þ þ t3 ðR_ T 0 SÞ Z ðR T 0 SÞ € GN theories predict similar patterns for the temperature and stresses. Moreover, it may be found from the figures b TÞ _ t3 r : ðK rTÞ r : ðC ¼ Z r : ðK rTÞ t3 r : ðK rTÞ _ that the GL theory predicts a Dirac-delta-type result for ð59Þ the stresses due to the existence of term T_ in the stress–strain relations. Also, the temperature patterns Also, the second term, which only contributes to the LS obtained from these theories are similar when a tempera- theory, should be substituted by using Eq. (6). Thus, for ture shock load is applied to the boundary of the sphere. simplicity, we only use the terms which appear in the LS theory. Eq. (6) for the LS theory, the case (Z ¼ 1; t3 ¼ 0), 5. Conclusions may be solved to find the heat flux in terms of temperature. When the initial heat flux through the body is zero we have In this paper, a new unified formulation for the Z t generalized coupled thermoelasticity theories based on q ¼ B 1 e s1 Kr T de B t. (60) the LS, GL, and GN models is proposed. The unifier 0 parameters are introduced to consolidate the equations of Differentiation of the above equation with respect to time the LS, GL, and GN theories in a single system of variable, considering that ðqðB1 Þ=qt ¼ B1 qB=qt B1 ¼ equations for the anisotropic and heterogeneous materials. s1 B1 , and multiplying the result by s leads to Then, the unified equations are reduced for the isotropic Z t and homogeneous materials. Transforming the governing s_q ¼ B1 e s1 Kr T de B t Kr T. (61) equations for the isotropic and homogeneous materials 0 ARTICLE IN PRESS A. Bagri, M.R. Eslami / International Journal of Mechanical Sciences 49 (2007) 1325–1335 1335 Substituting Eq. (61) in Eq. (59) gives [13] Nayfeh AH, Nemat-Nasser S. Thermoelastic waves in solids with Z t thermal relaxation. Acta Mechanica 1971;12:53–69. _ þ Z B1 Z ðR T 0 SÞ e s1 Kr T de B € t þ t3 ðR_ T 0 SÞ [14] Amin AM, Sierakowski RL. Effect of thermomechanical coupling of 0 the response of elastic solids. AIAA Journal 1990;28:1319–22. b TÞ. _ t3 r : ðK rTÞ r : ðC _ [15] Chen J, Dargush GF. Boundary element method for dynamic ¼ t3 r : ðK rTÞ ð62Þ poroelastic and thermoelastic analysis. International Journal of Now substituting S from Eq. (5) into Eq. (62) and Solids and Structures 1995;32(15):2257–78. [16] Chen H, Lin H. Study of transient coupled thermoelastic problems utilizing Eq. (2), considering the values for the unifier terms with relaxation times. Transactions of the ASME 1995;62:208–15. Z and t3 and the parameters t0 ; t1 ; t2 , and noting that b and [17] Hosseini Tehrani P, Eslami MR. Boundary element analysis of E are symmetric tensors, Eq. (7) is obtained. coupled thermoelasticity with relaxation times in finite domain. AIAA Journal 2000;38(3):534–41. References [18] Hosseini Tehrani P, Eslami MR. Boundary element analysis of finite domains under thermal and mechanical shock with the Lord- Shulman theory. Journal of Strain Analysis 2003;38(1):53–64. [1] Maxwell JC. On the dynamical theory of gases. Philosophical Transactions of the Royal Society of London 1867;157:49–88. [19] Bagri A, Eslami MR. Generalized coupled thermoelasticity of disks [2] Chester M. Second sound in solids. Physical Review 1963;131:2013–5. based on the Lord–Shulman model. Journal of Thermal Stresses [3] Lord HW, Shulman Y. A generalized dynamical theory of thermo- 2004;27(8):691–704. elasticity. Journal of the Mechanics and Physics of Solids [20] Chandrasekharaiah DS. A uniqueness theorem in the theory of 1967;15:299–309. thermoelasticity without energy dissipation. Journal of Thermal [4] Green AE, Lindsay KA. Thermoelasticity. Journal of Elasticity Stresses 1996;19:267–72. 1972;2(1):1–7. [21] Chandrasekharaiah DS. One-dimensional wave propagation in the [5] Green AE, Naghdi PM. A re-examination of the basic postulates of linear theory of thermoelasticity without energy dissipation. Journal thermomechanics. Proceedings of the Royal Society of London Series of Thermal Stresses 1996;19:695–710. A 1991;432:171–94. [22] Chandrasekharaiah DS. Complete solutions in the theory of [6] Green AE, Naghdi PM. Thermoelasticity without energy dissipation. thermoelasticity without energy dissipation. Mechanics Research Journal of Elasticity 1993;31:189–208. Communications 1997;24:625–30. [7] Ignaczak J. Linear dynamic thermoelasticity, a survey. Shock and [23] Sharma JN, Chauhan RS. Mechanical and thermal sources in a Vibration Digest 1981;13:3–8. generalized thermoelastic half-space. Journal of Thermal Stresses [8] Ignaczak J. Domain of influence results in generalized thermoelas- 2001;24:651–75. ticity, a survey. Applied Mechanics Reviews 1991;44(9):375–82. [24] Taheri H, Fariborz S, Eslami MR. Thermoelasticity solution of a [9] Francis PH. Thermomechanical effects in elastic wave propagation: a layer using the Green–Naghdi model. Journal of Thermal Stresses survey. Journal of Sound Vibration 1972;21:181–92. 2004;27(9):795–809. [10] Chandrasekharaiah DS. Thermoelasticity with second sound: a [25] Noda N, Furukawa T, Ashida F. Generalized theory in an infinite review. Applied Mechanics Reviews 1986;39(3):355–76. solid with a hole. Journal of Thermal Stresses 1989;12:385–402. [11] Chandrasekharaiah DS. Hyperbolic thermoelasticity: a review of [26] Abramovitz M, Stegun IA, editors. Handbook of mathematical recent literature. Applied Mechanics Reviews 1998;51(12, part functions, 9th Printing. New York: Dover; 1972, 1964. 1):705–29. [27] Honig G, Hirdes U. A method for the numerical inversion of Laplace [12] Nayfeh AH. Propagation of thermoelastic distribution in non- transforms. Journal of Computational and Applied Mathematics Fourier solids. AIAA Journal 1977;15:957–60. 1984;10:113–32.