Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 54, 1, pp. 295-304, Warsaw 2016 DOI: 10.15632/jtam-pl.54.1.295 APPLICATION OF FRACTIONAL ORDER THEORY OF THERMOELASTICITY TO A 1D PROBLEM FOR A SPHERICAL SHELL Waleed E. Raslan Mansoura University, Department of Mathematics and Engineering Physics, Mansoura, Egypt e-mail: w raslan@yahoo.com In this work, we apply the fractional order theory of thermoelasticity to a one-dimensional problem of distribution of thermal stresses and temperature in a generalized thermoelastic medium in the formof a spherical shell subjected to sudden change in the temperature of its external boundary. Laplace transform techniques are used to solve the problem. Numerical results are computed and represented graphically for the temperature, displacement and stress distributions. Keywords: fractional calculus, spherical shell, thermoelasticity 1. Introduction Biot (1956) formulated theory of coupled thermoelasticity to eliminate the paradox inherent in the classical uncoupled theory that elastic changes have no effect on temperature. Lord andShul- man (1967) introduced theory of generalized thermoelasticity with one relaxation time by using theMaxwell-Cattaneo law of heat conduction instead of the conventional Fourier 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. Sherief and El-Maghraby (2003, 2005) solved some crack problems for this theory. Sherief and Hamza (1994, 1996) obtained a solution to axisymmetric problems in spherical and cylindrical regions. Sherief and Ezzat (1994) obtained the solution in form of a series. Sherief et al. (2005) extended this theory to deal with micropolar materials. That theory was extended to deal with viscoelastic effects by Sherief et al. (2011). Lately, Sherief and Hussein (2012) developed theory of generalized poro-thermoelasticity. Fractional calculus has been successfully used to modify many existing models of physical processes, see Hilfer (2000), Sherief et al. (2012), Tenreiro et al. (2013). One can state that the whole theory of fractional derivatives and integrals was established in the 2nd half of the 19th century. A good review of the subject can be found in Podlubny (1998), Kaczorek (2011), Kaczorek and Rogowski (2015). Caputo and Mainardi (1971a,b) and Caputo (1974) found a good agreementwith experimental results bymaking use of fractional derivatives for description of viscoelastic materials and established the connection between the fractional derivatives and the theory of linear viscoelasticity. Adolfsson et al. (2005) constructed a new fractional order model of viscoelasticity. Povstenko (2009) made a review of thermoelasticity that uses a fractional heat conduction equation and proposed and investigated newmodels that incorporate fractional derivatives (Po- vstenko, 2005, 2011). Recently, the fractional order theory of thermoelasticity was derived by Sherief et al. (2010). It was a generalization of both coupled and generalized theories of thermo- elasticity. Some contributions to that theory are the works by Raslan (2015), Sherief and Abd El-Latief (2014a,b, 2015). 296 W.E. Raslan The aim of the present work is to solve a 1D problem for a spherical shell of a homogeneous, isotropic, thermoelastic medium occupying the region a ¬ r ¬ b subjected to thermal shock, using the fractional theory of thermoelasticity. The main reason behind the introduction of the fractional theory is that it predicts a retarded response to physical stimuli, as is found in nature and as opposed to the instantaneous response predicted by the generalized theory of thermoelasticity (Raslan, 2015). 2. Formulation of the problem In thiswork,we consider a 1Dproblem for a spherical shell of a homogeneous, isotropic, thermo- elastic medium occupying the region a¬ r¬ b, using the fractional theory of thermoelasticity. The outer surface of the shell is taken to be traction free and is subjected to thermal shock that is a function of time. The inner surface of the shell is thermally isolated by a rigid material. From physics of the problem, all functions will depend on the radial distance r and time t only. The displacement vector has only one non-zero component u(r,t) in the radial direction. Thegoverning equations, in the absenceof body forces andheat sources, are givenby (Sherief et al., 2010) (λ+2µ) ∂e ∂r −γ∂T ∂r = ρ ∂2u ∂t2 k∇2T = ( ∂ ∂t + τ0 ∂α+1 ∂tα+1 ) (ρcET +γT0e) σrr =λe+2µ ∂u ∂r −γ(T −T0) qr+ τ0 ∂αqr ∂tα =−k∂T ∂r (2.1) where T is the absolute temperature, ρ is density, λ and µ are Lamé’s constants and γ = αt(3λ+2µ), where αt is the coefficient of linear thermal expansion. T0 is the reference temperature assumed to be such that |(T − T0)/T0| ≪ 1 and α, τ0 are constants such that τ0 > 0, 0¬α¬ 1, cE is the specific heat per unit mass in the absence of deformation and k is the thermal conductivity, σrr is the normal stress component, qr is the component of the heat flux vector in the radial direction, and e is the cubical dilatation given by e= 1 r2 ∂ ∂r (r2u) (2.2) The operator ∇2 in the above equations is given by ∇2 = ∂2 ∂r2 + 2 r ∂ ∂r = 1 r2 ∂ ∂r ( r2 ∂ ∂r ) We shall use the following non-dimensional variables r∗ = cηr u∗ = cηu t∗ = c2ηt τ∗0 = c 2αηατ0 θ∗ = γ(T −T0) λ+2µ σ∗rr = σrr µ q∗ = γ k(λ+2µ) q where c= √ λ+2µ ρ η= ρcE k The governing equations, in non-dimensional form, are given by (with the asterisk dropped for convenience) ∂2u ∂t2 = ∂e ∂r − ∂θ ∂r ∇2θ= ( ∂ ∂t + τ0 ∂α+1 ∂tα+1 ) (θ+εe) σrr =(β 2−2)e+2∂u ∂r −β2θ qr+ τ ∂αqr ∂tα =−∂θ ∂r (2.3) Application of fractional order theory of thermoelasticity... 297 where ε= T0γ 2 λ+2µ kη β2 = λ+2µ µ In the above equation, the time fractional derivative of the order α used is taken to be in the sense of the Caputo fractional derivative. We assume that the boundary conditions have the form u(a,t) = 0 qr(a,t)= 0 σrr(b,t)= 0 θ(b,t)= f(t) (2.4) The initial conditions are taken to be homogeneous, i.e. we take u(r,t) ∣ ∣ ∣ t=0 = ∂u(r,t) ∂t ∣ ∣ ∣ ∣ ∣ t=0 =0 θ(r,t) ∣ ∣ ∣ t=0 = ∂θ(r,t) ∂t ∣ ∣ ∣ ∣ ∣ t=0 =0 σrr(r,t) ∣ ∣ ∣ t=0 = ∂σrr(r,t) ∂t ∣ ∣ ∣ ∣ ∣ t=0 =0 (2.5) 3. Solution in the Laplace transform domain Applying the Laplace transformwith the parameter s (denoted by the overbar) defined by the relation f(r,s)= ∞ ∫ 0 e−stf(r,t) dt (3.1) to both sides of equations (2.3), we get the following equations s2u= ∂e ∂r − ∂θ ∂r ∇2θ=(s+ τ0sα+1)(θ+εe) σrr = β2−2 r u+β2 ∂u ∂r −β2θ qr = −1 1+ τsα ∂θ ∂r (3.2) Applying the operator 1 r2 ∂ ∂r (r2 . . .) to equation (3.2)1, we obtain (∇2−s2)e=∇2θ (3.3) Eliminating θ between equations (3.2)2 and (3.3), we get { ∇4−∇2[s2+(1+ε)(s+ τ0sα+1)]+s3(1+τ0sα) } e=0 The above equation can be factorized as (∇2−k21)(∇2−k22)e=0 (3.4) where k21 and k 2 2 are the roots with positive real parts of the characteristic equation k4−k2[s2+(1+ε)(s+ τ0sα+1)]+s3(1+ τ0sα)= 0 (3.5) 298 W.E. Raslan where k21 and k 2 2 are given by k21 = s 2 { s+(1+ε)(1+ τ0s α)+ √ [s+(1+ε)(1+ τ0sα)]2−4s(1+ τ0sα) } k22 = s 2 { s+(1+ε)(1+ τ0s α)− √ [s+(1+ε)(1+ τ0sα)]2−4s(1+ τ0sα) } (3.6) Due to linearity, the solution to equation (3.4) has the form e= e1+e2 where ei is the solution to the following equation (∇2−k2i )ei =0 i=1,2 The above equation can be written as ∂2ei ∂r2 + 2 r ∂ei ∂r −k2iei =0 (3.7) Taking the substitution ei = gi√ r the above equation reduces to r2 ∂2gi ∂r2 +r ∂gi ∂r − ( k2r2+ 1 4 ) gi =0 This is the modified Bessel equation whose solution is gi =Aik 2 iI1/2(kir)+Bik 2 iK1/2(kir) Collecting the above results, the solution to (3.7) can be written as ei = 1√ r [Aik 2 iI1/2(kir)+Bik 2 iK1/2(kir)] (3.8) where Ai and Bi, i= 1,2 are parameters to be determined from the boundary conditions and Iµ(z), Kµ(z) are the modified Bessel functions of the first and second kinds of the order µ, respectively. Similarly, we can show that θi = 1√ r [A∗ik 2 iI1/2(kir)+B ∗ ik 2 iK1/2(kir)] (3.9) Substituting (3.8) and (3.9) into equation (3.3), we get A∗i =Ai(k 2 i −s2) B∗i =Bi(k2i −s2) (3.10) Substituting (3.10) into equation (3.9), one obtains θi = 1√ r [Ai(k 2 i −s2)I1/2(kir)+Bi(k2i −s2)K1/2(kir)] (3.11) Application of fractional order theory of thermoelasticity... 299 Thus we obtain e= 1√ r 2 ∑ i=1 [Aik 2 iI1/2(kir)+Bik 2 iK1/2(kir)] θ= 1√ r 2 ∑ i=1 [Ai(k 2 i −s2)I1/2(kir)+Bi(k2i −s2)K1/2(kir)] (3.12) Differentiating (3.12) with respect to r and substituting into (3.2)1, gives u= 1√ r 2 ∑ i=1 [AikiI3/2(kir)−BikiK3/2(kir)] (3.13) Differentiating (3.12)2 and (3.13) with respect to r, gives ∂θ ∂r = 1√ r 2 ∑ i=1 [Aiki(k 2 i −s2)I3/2(kir)−Biki(k2i −s2)K3/2(kir)] ∂u ∂r = 1√ r 2 ∑ i=1 { Aiki [ kiI1/2(kir)− 2 r I3/2(kir) ] +Biki [ kiK1/2(kir)+ 2 r K3/2(kir) ]} (3.14) Substituting (3.12) and (3.14)2 into equation (2.3)3, gives σrr = 1√ r 2 ∑ i=1 { Ai [ β2s2I1/2(kir)− 4 r kiI3/2(kir) ] +Bi [ β2s2K1/2(kir)+ 4 r kiK3/2(kir) ]} (3.15) Using equation (3.2)4, boundary conditions (2.4) can be written in the Laplace transform as u(a,s)= 0 ∂θ(a,s) ∂r =0 σrr(b,s)= 0 θ(b,s)= f(s) (3.16) Applying boundary conditions (3.16) into equations (3.12)2, (3.13), (3.14)1 and (3.15), gives 2 ∑ i=1 [AikiI3/2(kia)−BikiK3/2(kia)] = 0 2 ∑ i=1 [Aiki(k 2 i −s2)I3/2(kia)−Biki(k2i −s2)K3/2(kia)] = 0 2 ∑ i=1 { Ai [ β2s2I1/2(kib)− 4 b kiI3/2(kib) ] +Bi [ β2s2K1/2(kib)+ 4 b kiK3/2(kib) ]} =0 2 ∑ i=1 [Ai(k 2 i −s2)I1/2(kib)+Bi(k2i −s2)K1/2(kib)] = √ bf(s) The above equations can be put in the following form a11A1+a12B1+a13A2+a14B2 =0 a21A1+a22B1+a23A2+a24B2 =0 a31A1+a32B1+a33A2+a34B2 =0 a41A1+a42B1+a43A2+a44B2 = √ bf(s) 300 W.E. Raslan where a11 = k1I3/2(k1a) a12 =−k1K3/2(k1a) a13 = k2I3/2(k2a) a14 =−k2K3/2(k2a) a21 = k1(k 2 1 −s2)I3/2(k1a) a12 =−k1(k21 −s2)K3/2(k1a) a23 = k2(k 2 2 −s2)I3/2(k2a) a24 =−k2(k22 −s2)K3/2(k2a) a31 =β 2s2I1/2(k1b)− 4bk1I3/2(k1b) a32 =β 2s2K1/2(k1b)+ 4 b k1K3/2(k1b) a33 =β 2s2I1/2(k2b)− 4bk2I3/2(k2b) a34 =β 2s2K1/2(k2b)+ 4 b k2K3/2(k2b) a41 =(k 2 1 −s2)I1/2(k1b) a42 =(k21 −s2)K1/2(k1b) a43 =(k 2 2 −s2)I1/2(k2b) a44 =(k22 −s2)K1/2(k2b) Solving the above equations, we obtain A1 =− 1 Γ [a12(a23a34−a24a33)+a13(a24a32−a22a34)+a14(a22a33−a23a32)] √ bf(s) B1 = 1 Γ [a11(a23a34−a24a33)+a13(a24a31−a21a34)+a14(a21a33−a23a31)] √ bf(s) A2 =− 1 Γ [a11(a22a34−a24a32)+a12(a24a31−a21a34)+a14(a21a32−a22a31)] √ bf(s) B2 = 1 Γ [a11(a22a33−a23a32)+a12(a23a31−a21a33)+a13(a21a32−a22a31)] √ bf(s) where Γ = a11[a22(a33a44−a34a43)+a23(a34a42−a32a44)+a24(a32a43−a33a42)] −a12[a21(a33a44−a34a43)+a23(a34a41−a31a44)+a24(a31a43−a33a41)] +a13[a21(a32a44−a34a42)+a22(a34a41−a31a44)+a24(a31a42−a32a41)] −a14[a21(a32a43−a33a42)+a22(a33a41−a31a43)+a23(a31a42−a32a41)] 4. Inversion of the Laplace transform Let f(s) be the Laplace transform of f(t). The inversion formula for the Laplace transform has the form (Honig and Hirdes, 1984) f(t)= edt 2π ∞ ∫ −∞ eityf(d+iy) dy where d is a number greater than all the real parts of the singularities of f(s). Using Fourier series over the interval [0,2L], we get (Honig and Hirdes, 1984) f(t)≈ fN(t)= 1 2 c0+ N ∑ k=1 ck for 0¬ t¬ 2L (4.1) where ck = edt L Re [ e ikπt L f ( d+ ikπ L )] (4.2) The ‘Korrecktur’ method has been used to reduce the discretization error while the ε-algorithm has been used to reduce the truncation error (Honig and Hirdes, 1984). Application of fractional order theory of thermoelasticity... 301 5. Numerical results Copper has been chosen for purposes of numerical evaluations. The constants of the considered problem are shown in Table 1. Table 1 k=386W/(mK) αt =1.78 ·10−5K−1 cE =381J/(kgK) η=8886.73 µ=3.86 ·1010kg/(ms2) λ=7.76 ·1010kg/(ms2) ρ=8954kg/m3 T0 =293K ε=0.0168 τ0 =0.025s The computations have been carried out for a function f(t) given by f(t)=H(t) for which f(s)= 1 s The computations have been carried out for one value of time, namely t = 0.05, and two values of α, namely α=0.5 and α=1. The temperature, displacement and stress distributions have been obtained and plotted as shown in Figs. 1, 2 and 3, respectively. Fig. 1. Temperature distribution for t=0.05 Fig. 2. Displacement distribution for t=0.05 Fig. 3. Stress distribution for t=0.05 302 W.E. Raslan Next, the computations have been carried out for one value ofα, namely α=0.99, and two values of time, t = 0.05 and t = 0.1. The temperature, displacement and stress distributions have been obtained and plotted as shown in Figs. 4, 5 and 6, respectively. Fig. 4. Temperature distribution for α=0.99 Fig. 5. Displacement distribution for α=0.99 Fig. 6. Stress distribution for α=0.99 For the pervious steps, FORTRAN programming language has been used on a personal computer. Themaintained accuracy has been 5 digits for the numerical program. 6. Conclusions The computations show that: • For α=0.5, the solution behaves like in the coupled theory of thermoelasticity where the velocity of the wave is infinite, but forα=1 the solution becomes that of the generalized theory of thermoelasticity. • For α≈ 1, the solution seems to behave like in the generalized theory of thermoelasticity. This result is very important since the new theory may preserve the advantage of the generalized theory that the velocity of waves is finite. It is difficult to say whether the solution for α approaching 1 has a jump at the wave front or it is continuous with very fast changes (Povstenko, 2011). This aspect invites further investigation. Application of fractional order theory of thermoelasticity... 303 References 1. AdolfssonK., EnelundM., OlssonP., 2005,On the fractional ordermodel of viscoelasticity, Mechanics of Time Dependent Materials, 9, 15-34 2. Biot M., 1956, Thermoelasticity and irreversible thermo-dynamics, Journal of Applied Physics, 27, 240-253 3. Caputo M., 1974, Vibrations on an infinite viscoelastic layer with a dissipativememory, Journal of the Acoustical Society of America, 56, 897-904 4. Caputo M., Mainardi F., 1971a, A new dissipation model based onmemory mechanism, Pure and Applied Geophysics, 91, 134-147 5. Caputo M., Mainardi F., 1971a, Linear model of dissipation in anelastic solids, Rivis ta del Nuovo Cimento, 1, 161-198 6. Hilfer R., 2000,Applications of Fractional Calculus in Physics, World Scientific Publishing, Sin- gapore 7. Honig G., Hirdes U., 1984, A method for the numerical inversion of the Laplace transform, Journal of Computational and Applied Mathematics, 10, 113-132 8. Kaczorek T., 2011, Selected Problems of Fractional Systems Theory, Springer, Berlin 9. Kaczorek T., Rogowski K., 2015,Fractional Linear Systems and Electrical Circuits, Springer, Berlin 10. LordH., ShulmanY., 1967,A generalized dynamical theory of thermo-elasticity, Journal of the Mechanics and Physics of Solids, 15, 299-309 11. Podlubny I., 1998,Fractional Differential Equations, Academic Press, NY 12. Povstenko Y.Z., 2005, Fractional heat conduction and associated thermal stress, Journal of Thermal Stresses, 28, 83-102 13. Povstenko Y.Z., 2009, Thermoelasticity that uses fractional heat conduction equation, Journal of Mathematical Sciences, 162, 296-305 14. Povstenko Y.Z., 2011, Fractional Cattaneo-type equations and generalized thermoelasticity, Jo- urnal of Thermal Stresses, 34, 97-114 15. RaslanW., 2015,Application of fractional order theory of thermoelasticity in a thick plate under axisymmetric temperature distribution, Journal of Thermal Stresses, 38, 733-743 16. Sherief H., Abd El-Latief A.M., 2014a, Application of fractional order theory of thermoela- sticity to a 1D problem for a half-space, ZAMM, 94, 509-515 17. Sherief H., Abd El-Latief A.M., 2014b, Application of fractional order theory of thermoela- sticity to a 2D problem for a half-space,Applied Mathematics and Computation, 248, 584-592 18. Sherief H., Abd El-Latief A.M., 2015, A one-dimensional fractional order thermoelastic pro- blem for a spherical cavity,Mathematics and Mechanics of Solids, 20, 512-521 19. Sherief H., Allam M., El-Hagary M., 2011, Generalized theory of thermoviscoelasticity and a half-space problem, International Journal of Thermophysics, 32, 1271-1295 20. SheriefH.,El-MaghrabyN., 2003,An internal penny-shaped crack in an infinite thermoelastic solid, Journal of Thermal Stresses, 26, 333-352 21. SheriefH., El-MaghrabyN., 2005,Amode I crack problem for an infinite space in generalized thermoelasticity, Journal of Thermal Stresses, 28, 465-484 22. Sherief H., El-Sayed A.M.A., Abd El-Latief A.M., 2010, Fractional order theory of ther- moelasticity, International Journal of Solids and Structures, 47, 269-275 23. SheriefH., El-SayedA.M.A., Behiry S.H., RaslanW.E., 2012,Using fractional derivatives to generalize the Hodgkin-Huxleymodel, Fractional Dynamics and Control, Springer, 275-282 304 W.E. Raslan 24. Sherief H., Ezzat M., 1994, Solution of the generalized problem of thermoelasticity in the form of series of functions, Journal of Thermal Stresses, 17, 75-95 25. Sherief H., Hamza F., 1994, Generalized thermoelastic problem of a thick plate under axisym- metric temperature distribution, Journal of Thermal Stresses, 17, 435-452 26. Sherief H., Hamza F., 1996, Generalized two-dimensional thermoelastic problems in spherical regions under axisymmetric distributions, Journal of Thermal Stresses, 19, 55-76 27. Sherief H., Hamza F., El-Sayed A., 2005, Theory of generalized micropolar thermoelasticity and an axisymmetric half-space problem, Journal of Thermal Stresses, 28, 409-437 28. Sherief H., Hussein E., 2012, A mathematical model for short-time filtration in poroelastic media with thermal relaxation and two temperatures,Transport in Porous Media, 91, 199-223 29. Tenreiro J., Alexandra M., Trujillo J., 2013, Science metrics on fractional calculus deve- lopment since 1966,Fractional Calculus and Applied Analysis, 16, 479-500 Manuscript received March 22, 2014; accepted for print September 22, 2015