Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

3460 
 

https://internationalpubls.com 

Analysis of Three-Dimensional Non-Homogeneous Fractional order 

Thermoelastic Problem of Thick Rectangular Plate with Internal Heat 

Generation 

Yogesh U. Naner𝟏, Satish G. Khavale𝟐, Kishor. R. Gaikwad𝟑 

1Department of Mathematics, Shri Vitthal Rukhmini Arts, Commerce and Science College, Sawana, Yavatmal, 

Maharashtra, India-445205. 
2Department of Mathematics, Statistics and Computer Science, Chetana’s H. S. College of Commerce and Economics, 

SMT. K. C. College of Art, Bandra (East), Mumbai, Maharashtra, India - 400051. 
3Post Graduate Department of Mathematics, N. E. S., Science College, Nanded, Maharashtra, India - 431605. 

Email- ynaner@gmail.com, khavalesatish8@gmail.com , drkr.gaikwad@yahoo.in , satish.khavale@chetanacollege.in  

 

Article History: 

Received: 19-08-2025 

Revised: 24-09-2025 

Accepted: 15-10-2025 

Abstract:  This paper investigates the thermal behaviour of a three-dimensional thick 

rectangular plate governed by a fractional-order derivative. The plate, occupying the 

region                                    𝐷 = {(𝑥, 𝑦, 𝑧) ∈ 𝑅3: 0 ≤ 𝑥 ≤ 𝑎, 0 ≤ 𝑦 ≤ 𝑏, 0 ≤ 𝑧 ≤ 𝑐}, is 

initially at an arbitrary temperature distribution 𝑓(𝑥, 𝑦, 𝑧). For 𝑡 > 0, the plate 

experiences internal heat generation represented by 𝑔(𝑥, 𝑦, 𝑧, 𝑡) Btu/hr ft3, while all 

boundary surfaces are maintained at zero temperature. The governing model is 

formulated using the Caputo fractional derivative, and analytical solutions for 

temperature, displacement, and thermal stresses are derived through the integral 

transform technique. Numerical simulations are carried out using PTC Mathcad 

software, and the results are illustrated graphically. The study highlights the influence 

of the fractional-order parameter on heat conduction, displacement, and stress 

distribution, demonstrating the effectiveness of fractional-order thermoelastic models in 

capturing nonlocal thermal effects. 

Introduction: The study of thermoelastic behaviour in solid materials has gained 

significant importance due to its wide applications in aerospace, mechanical, and civil 

engineering. Classical thermoelastic theories often assume instantaneous heat 

propagation, which contradicts physical reality. To overcome this limitation, fractional-

order thermoelasticity has been introduced as an effective approach to model heat 

conduction with memory and nonlocal effects. In this paper, a three-dimensional non-

homogeneous thick rectangular plate with internal heat generation is analysed using the 

Caputo fractional derivative. The objective is to investigate how the fractional-order 

parameter (α) influences temperature, displacement, and stress distributions. The 

proposed model provides a more generalized framework that includes both the classical 

and wave-type heat conduction as special cases, thus offering a realistic representation 

of transient thermal behaviour in thick plates. 

Objectives:  

1. To investigate the thermal behaviour of a three-dimensional thick rectangular plate 

subjected to internal heat generation. 

2. To analyse the effect of the fractional-order parameter (α) on temperature, 

displacement, and thermal stresses. 

3. To demonstrate the usefulness of fractional-order thermoelastic models in capturing 

nonlocal and memory effects in heat conduction. 

Methods: The governing equations are derived using the Caputo fractional derivative 

for time-fractional heat conduction. Integral transform techniques are employed to 

obtain analytical solutions for temperature, displacement, and stress fields. Mittag-

mailto:ynaner@gmail.com
mailto:khavalesatish8@gmail.com
mailto:drkr.gaikwad@yahoo.in
mailto:satish.khavale@chetanacollege.in


Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

3461 
 

https://internationalpubls.com 

Leffler functions are used in the solutions to describe the fractional behaviour. 

Numerical simulations are performed in PTC Mathcad Prime for a thick rectangular 

plate using physical constants.  

Results: The temperature distribution decreases with increasing fractional order α and 

is maximum at the centre of the plate. Displacement components show symmetrical 

behaviour, increasing from the edges and vanishing at the centre. Stress components are 

compressive, peaking at the middle of the plate and zero at the edges. The fractional 

order parameter α significantly affects heat transfer. Graphs demonstrate clear variation 

between classical (α = 1) and fractional-order (α ≠ 1) thermoelastic behaviour.  

Conclusions: For α = 1, the model reduces to the classical heat diffusion equation, while 

α = 2 corresponds to the wave equation. The fractional order 0 < α < 1, α = 1, and 1 < α 

< 2 represent weak, normal, and strong conductivity, respectively. Stress distributions 

are tensile/compressive and follow normal curves in all directions, directly proportional 

to α. The fractional order α governs nonlocal heat transfer, influencing both response 

time and temperature overshoot. The model effectively describes transient thermoelastic 

behavior in thick plates and can be applied to other nonhomogeneous materials. 

Keywords: Thermal Behavior, Internal Heat Source, Caputo Fractional Derivative, 

Thick Rectangular Plate, Mittag-Leffler Functions. 

 

1. Introduction 

Biot [1] 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 are of the diffusion type predicting infinite speeds of propagation for heat waves 

contrary to physical observations. The generalised thermoelastic problem of a thick plate under axi-

symmetric heat supply was resolved by Sherief and Hamza [2]. The generalised dynamical theory of 

thermoelasticity with one relaxation period, for the isotropic body, was developed by Lord and 

Shulman [3]. In their study of the theoretical analysis of a rectangular plate in a compressive stress 

field and associated thermal stresses was resolved by Tanigawa and Komatsubara [4].   

The solution to the plane thermoelasticity problem for a rectangular domain is discussed in Vihak et 

al., [5]. A laminated rectangular plate subjected to a thermal shock was the focus of an investigation 

by Adam and Best [6]. Tanigawa et al. [7 8]. investigate the thermal analysis of thermal buckling 

caused by uniform heat supply in an orthotropic nonhomogeneous rectangular plate. The transient 

thermoelastic deformations of a thick functionally graded plate were investigated by Qian and Batra 

[9]. The thermal behaviour of thermoelastic thick plates under lateral loads was presented by Sharma 

et. al. [10]. Based on the equation of heat conduction in 1D and 2D with a time-fractional derivative 

and related thermal stresses, Povstenko [11] resolved a few thermoelastic issues. Some contribution of 

thermoelastic problems and fractional order thermoelastic problems have been discussed in [12-34]. 

2. Preliminaries 

In this section we give some definitions, notations and facts that we need in this article. 

 

• The Riemann-Liouville Fractional Differential Operator 

The definition of the Riemann-Liouville fractional differential operator given by [37]     



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

3462 
 

https://internationalpubls.com 

Suppose that 𝛼 > 0, t > a, 𝛼, a, t ∈ R. Then 

𝐷𝛼𝑓(𝑡) =

{
 
 

 
 1

𝛤(𝑛 − 𝛼)

𝑑𝑛

𝑑𝑡𝑛
∫

𝑓(𝜏)

(𝑡 − 𝜏)𝛼+1−𝑛
𝑑𝜏 

𝑡

0

, 𝑛 − 1 < 𝛼 < 𝑛 ∈ 𝑁

𝑑𝑛𝑓(𝑡)

𝑑𝑡𝑛
, 𝛼 = 𝑛 ∈ 𝑁

                                          (2.1) 

is called the Riemann-Liouville fractional derivative or the Riemann-Liouville fractional differential 

operator of order 𝛼. 

 

• The Caputo Fractional Differential Operator 

The definition of the Caputo fractional differential operator given by [37]  

Suppose that 𝛼 > 0, t > a, 𝛼, a, t ∈ R. Then 

𝐷𝛼𝑓(𝑡) =

{
 
 

 
 1

𝛤(𝑛 − 𝛼)

𝑑𝑛

𝑑𝑡𝑛
∫

𝑓(𝑛)(𝜏)

(𝑡 − 𝜏)𝛼+1−𝑛
𝑑𝜏 

𝑡

0

, 𝑛 − 1 < 𝛼 < 𝑛 ∈ 𝑁

𝑑𝑛𝑓(𝑡)

𝑑𝑡𝑛
, 𝛼 = 𝑛 ∈ 𝑁

                                          (2.2) 

is called the Caputo fractional derivative or the Caputo fractional differential operator of order 𝛼. 

 

• The Mittag-Leffler Function 

The Mittag-Leffler function is a generalization of the exponential function, first introduced as a one-

parameter function by the series  

𝐸𝛼(𝑧) = ∑
𝑧𝑘

𝛤(𝛼𝑘 + 1)

∞

𝑘=0

, 𝛼 > 0,  𝛼 𝜖 𝑅, 𝑧  ∈  𝐶                                                                                  (2.3) 

Later, the two-parameter generalization Mittag-Leffler function as: 

𝐸𝛼,𝛽(𝑧) =∑
𝑧𝑘

𝛤(𝛼𝑘 + 𝛽)

∞

𝑘=0

, 𝛼, 𝛽 > 0,  𝛼, 𝛽  ∈  𝑅, 𝑧  ∈  𝐶                                                                   (2.4) 

3. Problem formulation 

Consider a thick rectangular plate with length 𝑎, width 𝑏 and thickness 𝑐 occupying the space              

𝐷 = {(𝑥, 𝑦, 𝑧) ∈ 𝑅3: 0 ≤ 𝑥 ≤ 𝑎, 0 ≤ 𝑦 ≤ 𝑏, 0 ≤ 𝑧 ≤ 𝑐}. Initially the rectangular plate is at arbitrary 

temperature 𝑓(𝑥, 𝑦, 𝑧). For time 𝑡 > 0, heat is generated within the plate at a rate of 𝑔(𝑥, 𝑦, 𝑧, 𝑡) Btu/hr 

ft3, while the remaining boundaries are kept at zero temperature. An equation of the Caputo type time 

fractional differential equation of order 𝛼 is used to create a mathematical model that consider a thick 

rectangular plate with internal heat generation. 

The temperature 𝑇(𝑥, 𝑦, 𝑧, 𝑡) of the thick rectangular plate at time 𝑡 satisfying the time fractional 

heat conduction equation 

∂2𝑇

∂𝑥2
+
∂𝑇

∂𝑦2
+
∂2𝑇

∂𝑧2
+
𝑔(𝑥, 𝑦, 𝑧, 𝑡)

𝑘𝑡
=
1

𝑘

∂𝛼𝑇

∂𝑡𝛼
,                                                                                               (3.1) 

where 𝑘𝑡 and 𝑘 are thermal conductivity and thermal diffusivity of the material of the plate. 

with the boundary conditions, 

𝑇 = 0,  at all boundary surfaces for 𝑡 > 0,                                                                                          (3.2) 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

3463 
 

https://internationalpubls.com 

and initial conditions, 

𝑇 = 𝑓(𝑥, 𝑦, 𝑧),  at 𝑡 = 0, 0 < 𝛼 ≤ 1 ,                                                                                                 (3.3) 

∂𝑇

∂𝑡
= 0,  at 𝑡 = 0, 1 < 𝛼 ≤ 2                                                                                                   (3.4) 

 

 

Fig. 1. Geometrical representation of the plate. 

4. Analysis 

Here the plate is assumed sufficiently thick and considered free from traction. Since the plate is in a 

plane stress state without bending. Airy stress function method is applicable to the analytical 

development of the thermoelastic field. Airy stress function 𝑈(𝑥, 𝑦, 𝑧, 𝑡) which satisfy the following 

relation 

(
∂2

∂𝑥2
+

∂

∂𝑦2
+
∂2

∂𝑧2
)

2

𝑈 = −𝜆𝐸 (
∂2

∂𝑥2
+

∂

∂𝑦2
+
∂2

∂𝑧2
)𝑇                                                                         (4.1) 

where 𝜆 and 𝐸 are linear coefficient of the thermal expansion, Youngs modulus elasticity of the 

material of the plate. 

The displacement components 𝑢𝑥, 𝑢𝑦 and 𝑢𝑧 in the 𝑋, 𝑌 and 𝑍 direction are represented in the 

integral form as 

𝑢𝑥 = ∫[
1

𝐸
(
∂2𝑈

∂𝑦2
+
∂2𝑈

∂𝑧2
− 𝜈

∂2𝑈

∂𝑥2
) + 𝜆𝑇] 𝑑𝑥                                                                                           (4.2) 

𝑢𝑦 = ∫[
1

𝐸
(
∂2𝑈

∂𝑥2
+
∂2𝑈

∂𝑧2
− 𝜈

∂2𝑈

∂𝑦2
) + 𝜆𝑇] 𝑑𝑦                                                                                            (4.3) 

𝑢𝑧 = ∫[
1

𝐸
(
∂2𝑈

∂𝑥2
+
∂2𝑈

∂𝑦2
− 𝜈

∂2𝑈

∂𝑧2
) + 𝜆𝑇] 𝑑𝑧                                                                                            (4.4) 

where 𝜈 is the poissons ratio of the material of the plate. 

The stress components are: 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

3464 
 

https://internationalpubls.com 

𝜎𝑥𝑥 = (
∂2𝑈

∂𝑦2
+
∂2𝑈

∂𝑧2
)                                                                                                                                    (4.5) 

𝜎𝑦𝑦 = (
∂2𝑈

∂𝑥2
+
∂2𝑈

∂𝑧2
)                                                                                                                                      (4.6) 

𝜎𝑧𝑧 = (
∂2𝑈

∂𝑥2
+
∂2𝑈

∂𝑦2
)                                                                                                                                     (4.7) 

Equations [3.1] to [4.7] constitute the mathematical formulation of the problem under consideration. 

 

5. Determination of Temperature Field 

To obtain the expression for temperature function 𝑇(𝑥, 𝑦, 𝑧, 𝑡); we introduce the ‘‘triple integral transform’’and 

its corresponding ‘triple-inversion formula’ as defined in [35] respectively as 

𝑇(𝛽𝑚, 𝜈𝑛, 𝜂𝑝, 𝑡) = ∫ ∫ ∫ 𝐾
𝑐

𝑧′=0

𝑏

𝑦′=0

𝑎

𝑥′=0

(𝛽𝑚, 𝑥
′). 𝐾(𝜈𝑛, 𝑦

′). 𝐾(𝜂𝑝, 𝑧
′)                                                  (5.1)

× 𝑇(𝑥′, 𝑦′, 𝑧′, 𝑡)𝑑𝑥′𝑑𝑦′𝑑𝑧′
 

𝑇(𝑥, 𝑦, 𝑧, 𝑡) = ∑∑∑𝐾

∞

𝑝=0

∞

𝑛=0

∞

𝑚=0

(𝛽𝑚, 𝑥). 𝐾(𝜈𝑛, 𝑦).𝐾(𝜂𝑝, 𝑧)𝑇(𝛽𝑚, 𝜈𝑛, 𝜂𝑝 , 𝑡)                                        (5.2)
 

where the kernels 

𝐾(𝛽𝑚, 𝑥) = √
2

𝑎
sin(𝛽𝑚𝑥)                                                                                                                              (5.3) 

𝐾(𝜈𝑛, 𝑦) = √
2

𝑏
sin(𝜈𝑛𝑦)                                                                                                                                (5.4) 

𝐾(𝜂𝑝, 𝑧) = √
2

𝑐
sin(𝜂𝑝𝑧)                                                                                                                                (5.5) 

and eigenvalues are 

𝛽𝑚 is 𝑚𝑡ℎ root of transcendental equation 

sin(𝛽𝑚 . 𝑎) = 0 

i.e. 

𝛽𝑚 =
𝑚𝜋

𝑎
,  𝑚 = 1,2,3,…                                                                                                                      (5.6) 

𝜈𝑛 is 𝑛𝑡ℎ root of transcendental equation 

sin(𝜈𝑛. 𝑏) = 0 

i.e. 

𝜈𝑛 =
𝑛𝜋

𝑏
,  𝑛 = 1,2,3,…                                                                                                                         (5.7) 

𝜂𝑝 is 𝑝𝑡ℎ root of transcendental equation 

sin(𝜂𝑝. 𝑐) = 0 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

3465 
 

https://internationalpubls.com 

i.e. 

𝜂𝑝 =
𝑝𝜋

𝑐
,  𝑝 = 1,2,3, …                                                                                                                       (5.7) 

On applying triple-integral transform defined in equation [5.1] to equations [3.1]–[3.4] and then 

using their inversions defined in equation [4.2], one obtains the expressions of the temperature as: 

𝑇(𝑥, 𝑦, 𝑧, 𝑡) = ∑∑∑𝐾

∞

𝑝=0

∞

𝑛=0

∞

𝑚=0

(𝛽𝑚, 𝑥). 𝐾(𝜈𝑛, 𝑦).𝐾(𝜂𝑝, 𝑧)
1

𝑘(𝛽𝑚2 + 𝜈𝑛2 + 𝜂𝑝2)

× [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)] [𝑓(𝛽𝑚, 𝜈𝑛, 𝜂𝑝) +

𝑘

𝑘𝑡

×∫ [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)]

𝑡

𝑡′=0

𝑔(𝛽𝑚 , 𝜈𝑛, 𝜂𝑝, 𝑡
′)]

                                     (5.9) 

where 𝐸𝛼,𝛼(. ) - two-parameter Mittag-Leffler function, 

𝑓(𝛽𝑚, 𝜈𝑛, 𝜂𝑝) = ∫ ∫ ∫ sin
𝑐

𝑧′=0

𝑏

𝑦′=0

𝑎

𝑥′=0

(𝛽𝑚𝑥
′)sin(𝜈𝑛𝑦

′)sin(𝜂𝑝𝑧
′)𝑓(𝑥′, 𝑦′, 𝑧′)𝑑𝑥′𝑑𝑦′𝑑𝑧′

𝑔(𝛽𝑚, 𝜈𝑛, 𝜂𝑝, 𝑡
′) = ∫ ∫ ∫ sin

𝑐

𝑧′=0

𝑏

𝑦′=0

𝑎

𝑥′=0

(𝛽𝑚𝑥
′)sin(𝜈𝑛𝑦

′)sin(𝜂𝑝𝑧
′)𝑔(𝑥′, 𝑦′, 𝑧′, 𝑡′)𝑑𝑥′𝑑𝑦′𝑑𝑧′

 

6. Distributions of Displacement and Stresses 

Using equation [5.9] in [4.1], one obtains 

𝑈 =
8𝐸

𝑎𝑏𝑐
∑ ∑∑sin

∞

𝑝=0

∞

𝑛=0

∞

𝑚=0

(𝛽𝑚𝑥)sin(𝜈𝑛𝑦)sin(𝜂𝑝𝑧)
1

𝑘(𝛽𝑚2 + 𝜈𝑛2 + 𝜂𝑝2)
2

× [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)] [𝑓(𝛽𝑚, 𝜈𝑛, 𝜂𝑝) +

𝑘

𝑘𝑡

×∫ [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)]

𝑡

𝑡′=0

𝑔(𝛽𝑚, 𝜈𝑛, 𝜂𝑝 , 𝑡
′)]

                                           (6.1) 

Now using equations [5.9] and [6.1] in equations [4.2] to [4.4], one obtains the expressions for 

displacement as 

𝑢𝑥 =
−8

𝑎𝑏𝑐
∑ ∑∑(

cos(𝛽𝑚𝑥)

𝛽𝑚
)

∞

𝑝=0

∞

𝑛=0

∞

𝑚=0

sin(𝜈𝑛𝑦)sin(𝜂𝑝𝑧) [
(𝜈 − 1)𝛽𝑚

2 − 2(𝜈𝑛
2 + 𝜂𝑝

2)

𝑘(𝛽𝑚2 + 𝜈𝑛2 + 𝜂𝑝2)
2 ]

× [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)] [𝑓(𝛽𝑚, 𝜈𝑛, 𝜂𝑝) +

𝑘

𝑘𝑡

×∫ [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)]

𝑡

𝑡′=0

𝑔(𝛽𝑚 , 𝜈𝑛, 𝜂𝑝, 𝑡
′)]

                        (6.2) 

𝑢𝑦 =
−8

𝑎𝑏𝑐
∑ ∑∑sin

∞

𝑝=0

∞

𝑛=0

∞

𝑚=0

(𝛽𝑚𝑥)(
cos(𝜈𝑛𝑦)

𝜈𝑛
)sin(𝜂𝑝𝑧) [

(𝜈 − 1)𝜈𝑛
2 − 2(𝛽𝑚

2 + 𝜂𝑝
2)

𝑘(𝛽𝑚2 + 𝜈𝑛2 + 𝜂𝑝2)
2 ]

× [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)] [𝑓(𝛽𝑚, 𝜈𝑛, 𝜂𝑝) +

𝑘

𝑘𝑡

× ∫ [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)]

𝑡

𝑡′=0

𝑔(𝛽𝑚, 𝜈𝑛, 𝜂𝑝, 𝑡
′)]

                      (6.3) 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

3466 
 

https://internationalpubls.com 

𝑢𝑧 =
−8

𝑎𝑏𝑐
∑ ∑∑sin

∞

𝑝=0

∞

𝑛=0

∞

𝑚=0

(𝛽𝑚𝑥)sin(𝜈𝑛𝑦)(
cos(𝜂𝑝𝑧)

𝜂𝑝
) [
(𝜈 − 1)𝜂𝑝

2 − 2(𝛽𝑚
2 + 𝜈𝑛

2)

𝑘(𝛽𝑚2 + 𝜈𝑛2 + 𝜂𝑝2)
2 ]

× [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)] [𝑓(𝛽𝑚, 𝜈𝑛, 𝜂𝑝) +

𝑘

𝑘𝑡

×∫ [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)]

𝑡

𝑡′=0

𝑔(𝛽𝑚 , 𝜈𝑛, 𝜂𝑝, 𝑡
′)]

                        (6.4) 

Now using equations [5.9] and [6.1] in equations [4.5] to [4.7], one obtains expressions for thermal 

stresses as 

𝜎𝑥𝑥 =
−8𝐸

𝑎𝑏𝑐
∑ ∑∑sin

∞

𝑝=0

∞

𝑛=0

∞

𝑚=0

(𝛽𝑚𝑥)sin(𝜈𝑛𝑦)sin(𝜂𝑝𝑧)
(𝜈𝑛

2 + 𝜂𝑝
2)

𝑘(𝛽𝑚2 + 𝜈𝑛2 + 𝜂𝑝2)
2

× [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)] [𝑓(𝛽𝑚, 𝜈𝑛, 𝜂𝑝) +

𝑘

𝑘𝑡

× ∫ [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)]

𝑡

𝑡′=0

𝑔(𝛽𝑚, 𝜈𝑛, 𝜂𝑝, 𝑡
′)]

                                          (6.5) 

𝜎𝑦𝑦 =
−8𝐸

𝑎𝑏𝑐
∑ ∑∑sin

∞

𝑝=0

∞

𝑛=0

∞

𝑚=0

(𝛽𝑚𝑥)sin(𝜈𝑛𝑦)sin(𝜂𝑝𝑧)
(𝛽𝑚

2 + 𝜂𝑝
2)

𝑘(𝛽𝑚2 + 𝜈𝑛2 + 𝜂𝑝2)
2

× [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)] [𝑓(𝛽𝑚, 𝜈𝑛, 𝜂𝑝) +

𝑘

𝑘𝑡

×∫ [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)]

𝑡

𝑡′=0

𝑔(𝛽𝑚, 𝜈𝑛, 𝜂𝑝, 𝑡
′)]

                                           (6.6) 

𝜎𝑧𝑧 =
−8𝐸

𝑎𝑏𝑐
∑ ∑∑sin

∞

𝑝=0

∞

𝑛=0

∞

𝑚=0

(𝛽𝑚𝑥)sin(𝜈𝑛𝑦)sin(𝜂𝑝𝑧)
(𝛽𝑚

2 + 𝜈𝑛
2)

𝑘(𝛽𝑚2 + 𝜈𝑛2 + 𝜂𝑝2)
2

× [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)] [𝑓(𝛽𝑚, 𝜈𝑛, 𝜂𝑝) +

𝑘

𝑘𝑡

×∫ [𝑡𝛼−1𝐸𝛼,𝛼(−𝑘(𝛽𝑚
2 + 𝜈𝑛

2 + 𝜂𝑝
2)𝑡𝛼)]

𝑡

𝑡′=0

𝑔(𝛽𝑚, 𝜈𝑛, 𝜂𝑝, 𝑡
′)]

                                            (6.7) 

 

7. Numerical Calculations and Discussion 

Setting 

𝑓(𝑥, 𝑦, 𝑧) = 𝑒𝑡(𝑎 − 𝑥)(1 − 𝑒𝑥)(𝑏 − 𝑦)(1 − 𝑒𝑦)(𝑐 − 𝑧)(1 − 𝑒𝑧), 

𝑔(𝑥, 𝑦, 𝑧, 𝑡) = 𝑔𝑖𝛿(𝑥 − 𝑥1)𝛿(𝑦 − 𝑦1)𝛿(𝑧 − 𝑧1)𝛿(𝑡 − 𝜏) where 𝛿 is the Derac-delta function and 

𝑔(𝑥, 𝑦, 𝑧, 𝑡) is an instantaneous line heat source situated at the center of the rectangular plates. 

Dimension 

Length of rectangular plate 𝑎 = 3 ft 

Breadth of rectangular plate 𝑏 = 2 ft 

Height of rectangular plate 𝑐 = 1 ft 

Central length of rectangular plate 𝑥1 =1.5 ft 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

3467 
 

https://internationalpubls.com 

Central breadth of rectangular plate 𝑦1 = 1 ft 

Central 1eight of rectangular plate 𝑧1 = 0.5 ft 

For the purpose of doing a numerical calculation for a thick rectangular plate, the material aluminium 

(Pure) was selected. PTC Mathcad Prime, a computational mathematical software, was used to create 

the numerical calculations and graphics. 

The aluminium (pure) material was chosen for purposes of numerical evaluations. The parameters of 

the problem are given in FPS units, as illustrated in Table 1. 

TABLE -1:-  MATERIAL CONSTANT 

Physical constant Value 

Thermal diffusivity (k) 3.33 ft2/hr 

Thermal conductivity (𝑘𝑡) 117 Btu/(hr.ft.0F) 

Density (𝜌) 169 lb/ft3 

Specific heat (𝑐𝑝) 0.208 Btu/(lb.0F) 

Coefficient of linear thermal expansion (𝑎𝑡) 12.84×10−6 1/F 

Poisson ratio (𝜈) 0.35 

Lame’s constants (𝜇) 26.67 

Young’s modulus (𝐸) 70 GPa 

 

In this section, we analysed a fractional heat conduction problem for a thick rectangular plate. We 

have discussed the variation of temperature, displacement and thermal stresses for different values of 

fractional-order parameter 𝛼 = 0.5,1,1.5,2. The graphs are plotted for different fractional-order 

parameter 𝛼 = 0.5,1,1.5,2 at time 𝑡 = 1 hr. 

Figure 2 shows the temperature distribution for different values of fractional-order parameter 𝛼 =

0.5,1,1.5,2 in 𝑋-direction. It is observed that, the magnitude of the temperature increases with decrease 

the value of fractional-order parameter 𝛼 and it attains its peak at 𝑥 = 1.5. Also, it is zero at both edges 

(𝑥 = 0 and 𝑥 = 3) of the rectangular plate along the 𝑋-direction and it shows the normal curve. 

Figures 3–5 show the displacement distribution functions 𝑢𝑥, 𝑢𝑦 and 𝑢𝑧 in 𝑋, 𝑌 and 𝑍-directions for 

different values of 𝛼. It is clear that, the displacement functions increase from initial edge towards the 

extreme edge and it becomes zero at the middle part for different values of 𝛼. 

Figure 6-8 shows the stress distribution 𝜎𝑥𝑥, 𝜎𝑦𝑦 and 𝜎𝑧𝑧 for different values of 𝛼 in 𝑋, 𝑌 and 𝑍-

direction respectively. We observe that, the magnitude of both the stresses is maximum at the middle 

part with increasing the value of fractional-order parameter 𝛼 and becomes zero at the initial and 

extreme edges in the 𝑋, 𝑌 and 𝑍-direction respectively. Also, the stress components 𝜎𝑥𝑥, 𝜎𝑦𝑦 and 𝜎𝑧𝑧 

are compressive throughout the plate and shows the normal curve. 

 

 

  



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

3468 
 

https://internationalpubls.com 

 
Fig. 2. Temperature distribution for different values of 𝛼 in 𝑋-direction. 

 

 
Fig. 3. Displacement distribution 𝑢𝑥 for different values of 𝛼 in 𝑋-direction. 

 

 
 

Fig. 4. Displacement distribution 𝑢𝑦 for different values of 𝛼 in 𝑌-direction. 

 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

3469 
 

https://internationalpubls.com 

 
Fig. 5. Displacement distribution 𝑢𝑧 for different values of 𝛼 in 𝑍-direction. 

 

 

 
Fig. 6. Stress distribution 𝜎𝑥𝑥 for different values of 𝛼 in 𝑋-direction. 

 

 

 
Fig. 7. Stress distribution 𝜎𝑦𝑦 for different values of 𝛼 in 𝑌-direction. 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

3470 
 

https://internationalpubls.com 

 

 
Fig. 8. Stress distribution 𝜎𝑧𝑧 for different values of 𝛼 in 𝑍-direction. 

 

8. Conclusion 

In this article, we proposed the time-fractional heat conduction equation with zero initial conditions. 

Integral transform technique is adopted for computing the temperature of heat conduction equation. 

The results are obtained in the form of two-parameter Mittag-Leffler function. Figures 2-8 shows the 

behavior of temperature, displacement and thermal stresses along 𝑋, 𝑌 and 𝑍 directions for the different 

values of the fractional-order parameter 𝛼 = 0.5,1,1.5,2 and shows the variation between classical and 

fractional-order thermoelasticity. 

The following concluding remakes can be considered according to the results of the present study. 

1. When 𝛼 = 1, all the graphical representation satisfied the heat diffusion equation, whereas 𝛼 =

2, gives the representation of wave equations. 

2. It is observed that, for the different values of the fractional-order parameter 𝛼 in the 𝑋, 𝑌 and 

𝑍 directions the effect temperature, displacement and thermal stresses represents the weak, 

normal and strong conductivity, within the range of 0 < 𝛼 < 1, 𝛼 = 1 and 1 < 𝛼 < 2 

respectively. 

3. The stress distribution functions 𝜎𝑥𝑥, 𝜎𝑦𝑦 and 𝜎𝑧𝑧 are tensile and they form normal curve along 

the 𝑋, 𝑌 and 𝑍 direction respectively. Also, they are directly proportional to the different values 

of 𝛼 in the 𝑋, 𝑌 and 𝑍 direction respectively. 

4. Fractional order 𝛼 governs the non-local behavior of heat transfer process, and then alters the 

response time and overshooting occurrence of the temperature. Larger the parameter 𝛼, the 

slower the velocity of heat waves-like. Its effect is limited to short times and disappears for 

long time. 

Refrences 

[1] M. A. Biot., “Thermoelasticity and irreversible thermodynamics”, J. Appl. Phys, vol. 27, pp. 240–

253, 1956. 

[2] H. H. Sherief, F. A. Hamza, ‘‘Generalized thermoelastic problem of a thick plate under axi-

symmetric temperature distribution’’, Journal of Thermal Stresses, vol. 17(3), pp. 435–452, 1994. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

3471 
 

https://internationalpubls.com 

[3] H. W. Lord, Y. Shulman, ‘‘A generalized dynamical theory of thermoelasticity’’, J. Mech. Phys. 

Solids, vol. 15, pp. 299–307, 1967. 

[4] Y. Tanigawa, Y. Komatsubara, ‘‘Theoretical analysis of a rectangular plate and its thermal stress 

intensity factor compressive stress field’’, Journal of Thermal Stresses, vol. 20, pp. 517–542, 

1997. 

[5] V. M. Vihak, M. Y. Yuzvyak, A. V. Yasinskij, ‘‘The solution of the plane thermoelasticity 

problem for a rectangular domain’’, Journal of Thermal Stresses, vol. 21, pp. 545–561, 1998. 

[6] R. J.Adam, C. W. Best, ‘‘Thermoelastic vibrations of a laminated rectangular plate subjected to a 

thermal shock’’, Journal of Thermal Stresses, vol. 22, pp. 875–895, 1999. 

[7] T. Morimoto, Y. Tanigawa, R. Kawamura, ‘‘Thermal buckling analysis of an orthotropic 

nonhomogeneous rectangular plate due to uniform heat supply’’, Nippon Kikai Gakkai Zairyo 

Rikigaku Bumon Koenkai Koen Ronbunshu, vol. 202, pp. 283–284, 2002. 

[8] Y. Ootao, Y. Tanigawa, ‘‘Transient thermoelastic problem of functionally graded thick strip due 

to nonuniform heat supply’’, Composite Structures, vol. 63, pp. 139–146, 2004. 

[9] L. F. Qian, R. C. Batra, ‘‘Transient thermoelastic deformations of a thick functionally graded 

plate’’, Journal of Thermal Stresses, vol. 27, pp. 705–740, 2004. 

[10] J. N. Sharma , P. K. Sharma, R. L. Sharma, ‘‘Behavior of Thermoelastic Thick Plate under Lateral 

Loads’’, Journal of Thermal Stresses, vol. 27, pp. 171–191, 2004. 

[11] Y. Z. Povstenko, ‘‘Fractional heat conduction equation and associated thermal stresses’’, Journal 

of Thermal Stresses, 28, 83–102, 2005. 

[12] K. R. Gaikwad, K. P. Ghadle, ‘‘Quasi-static thermal stresses in a thick rectangular plate’’, Global 

Journal of Pure and Applied Mathematics vol. 5(2), 109-117, 2009. 

[13] H. H. Sherief, A. El-Said, A. Abd El-Latief, ‘‘Fractional order theory of thermoelasticity’’, 

International Journal of Solids and Structures, vol. 47, pp. 269–275, 2010. 

[14] K. R. Gaikwad, K. P. Ghadle, ‘‘Three dimensional non-homogeneous thermoelastic problem in a 

thick rectangular plate due to internal heat generation’’, Southern Africa Journal of Pure and 

Applied Mathematics, vol. 5, pp. 26–38, 2011. 

[15] K. R. Gaikwad, K. P. Ghadle, ‘‘Non-homogeneous heat conduction problem and its thermal 

deflection due to internal heat generation in a thin hollow circular disk’’, Journal of Thermal 

stresses, vol. 35, pp. 485–498, 2012. 

[16] K. R. Gaikwad, ‘‘Analysis of thermoelastic deformation of a thin hollow circular disk due to 

partially distributed heat supply’’, Journal of Thermal stresses, vol. 36, pp. 207–224, 2013. 

[17] K. R. Gaikwad, ‘‘Mathematical modelling and its simulation of a quasi-static thermoelastic 

problem in a semi-infinite hollow circular disk due to internal heat generation’’, Journal of the 

Korean Society for Industrial and Applied Mathematics, vol. 19, pp. 69–81, 2015. 

[18] V. R. Manthena, N. K. Lamba, G. D. Kedar, K. C. Deshmukh, ‘‘Effects of stress resultants on 

thermal stresses in a functionally graded rectangular plate due to temperature dependent material 

properties’’, Int. J. Thermodyn., vol. 19(4), pp. 235-242, 2016. 

[19] V. R. Manthena, N. K. Lamba, G. D. Kedar, ‘‘Transient thermoelastic problem of a 

nonhomogeneous rectangular plate’’, Journal Thermal Stresses, vol. 40(5), pp. 627-640, 2017. 

[20] K. R. Gaikwad, ‘‘Axi-symmetric thermoelastic stress analysis of a thin circular plate due to heat 

generation’’, International Journal of Dynamical Systems and Differential Equations, vol. 9, pp. 

187–202, 2019. 

[21] K. R. Gaikwad, S. G. Khavale, ‘‘Time fractional heat conduction problem in a thin hollow circular 

disk and it’s thermal deflection’’, Easy Chair, 1672, 1–11, 2019. 

[22] F. Yekkalam Tash, B. Navayi Neya, ‘‘An analytical solution for bending of transversely isotropic 

thick rectangular plates with variable thickness’’, Appl. Math. Model. vol. 77(2), pp. 1582-1602, 

2020. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 10s (2025) 

  

3472 
 

https://internationalpubls.com 

[23] S. G. Khavale, K. R. Gaikwad, ‘‘Generalized theory of magneto-thermo-viscoelastic Spherical 

cavity problem under Fractional order derivative: State Space Approach’’, Advances in 

Mathematics: Scientific Journal, vol. 9, pp. 9769–9780, 2020. 

[24] K. R. Gaikwad, S. G. Khavale, ‘‘Time fractional 2D thermoelastic problem of thin hollow circular 

disk and it’s associated thermal stresses’’, Bulletin of the Marathwada Mathematical Society, vol. 

21(1 & 2), pp. 37-47, 2020. 

[25] F. Y. Tash, B. N. Neya, ‘‘An analytical solution for bending of transversely isotropic thick 

rectangular plates with variable thickness’’, Appl. Math. Model., vol. 77(2), pp. 1582-– 1602, 

2020. 

[26] K. R. Gaikwad, Y. U. Naner, ‘‘Analysis of transient thermoelastic temperature distribution of a 

thin circular plate and its thermal deflection under uniform heat generation’’, Journal of Thermal 

Stress, vol. 44(1), pp. 75–85, 2021. 

[27] K. R. Gaikwad, Y. U. Naner, S. G. Khavale, ‘‘Time fractional thermoelastic stress analysis of a 

thin rectangular plate’’, NOVYI MIR Research Journal, vol. 6(1), pp. 42–56, 2021. 

[28] S. G. Khavale, K. R. Gaikwad, ‘‘Fractional order thermoelastic problem of thin hollow circular 

disk and its thermal stresses under axi-symmetric heat supply’’, Design Engineering, vol. 2021(9), 

pp. 13851–13862, 2021. 

[29] K. R. Gaikwad, V. G. Bhandwalkar, ‘‘Fractional order thermoelastic problem for finite 

piezoelectric rod subjected to different types of thermal loading - direct approach’’, Journal of the 

Korean Society for Industrial and Apllied Mathematics, vol. 25(3), pp. 117-131, 2021. 

[30] S. G. Khavale, K. R. Gaikwad, ‘‘Analysis of non-integer order thermoelastic temperature 

distribution and thermal deflection of thin hollow circular disk under the axi-symmetric heat 

supply’’, Journal of the Korean Society for Industrial and Applied Mathematics, vol. 26(1), pp. 

67-75, 2022. 

[31] K. R. Gaikwad, S. G. Khavale, ‘‘Fractional order transient thermoelastic stress analysis of a thin 

circular sector disk’’, International Journal of Thermodynamics, vol. 25(1), pp. 1-8, 2022. 

[32] S. G. Khavale, K. R. Gaikwad, ‘‘2D problem for a sphere in the fractional order theory 

thermoelasticity to axisymmetric temperature distribution’’, Advances in Mathematics: Scientific 

Journal, vol. 11(1), pp. 1–15, 2022. 

[33] S. G. Khavale, K. R. Gaikwad, ‘‘Two-dimensional generalized magneto-thermo-viscoelasticity 

problem for a spherical cavity with one relaxation time using fractional derivative’’, International 

Journal of Thermodynamics, vol. 25(2), pp. 89-97, 2022. 

[34] K. R. Gaikwad, Y. U. Naner, S. G. Khavale, ‘‘Transient thermoelastic bending analysis of a 

rectangular plate with a simply supported edge under heat source: Green’s function approach’’, 

Int. J. Nonlinear Anal. Appl., In Press, pp. 1-15, 2022. 

[35] N. M. Ozisik, Boundary Value Problem of Heat Conduction, International Textbook Company, 

Scranton, Pennsylvania, 1968. 

[36] I. N. Sneddon, The use of Integral Transform, McGraw Hill, New York, 1972. 

[37] I. Podlubny, Fractional differential Equation, Academic Press, San Diego, 1999. 

[38] B. G. Korenev, Bessel Functions and Their Applications. Boca Raton: CRC Press, 2003. 

[39] Y. Z. Povstenko, Fractional Thermoelasticity, New York: Springer, 2015. 

[40] PTC Mathcad Prime-7.0.0.0, [Online]. Available: https://support.ptc.com/help/mathcad/r7.0/en/ 

(accessed Sep. 1, 2022).