Fractional order model of micropolar thermoelasticity and 2D half-space problem

In this manuscript, the generalized micropolar theory of thermoelasticity is modified using fractional calculus. The revised equations are used to solve a problem for a half-space whose boundary is rigidly fixed and subjected to an axisymmetric thermal shock. Laplace and Hankel transform techniques are used. The analytical solution in the transform domain is obtained by using a new direct approach without the customary use of potential functions. By using a numerical method based on the Fourier expansion technique, the inverse of the double transform can be obtained. The numerical results for displacement, microrotation, stress, micro-stress, and temperature are obtained and represented graphically. Comparisons are made with the results of the older theory.


List of symbols e i j
Components of the strain tensor e e kk Cubical dilatation σ i j Components of the stress tensor t Time variable T Absolute temperature T 0 The temperature of the medium in its natural state, assumed to be such at T −T 0 Mass density c E Specific heat at constant strain k, k i j Thermal conductivity α t Coefficient of linear thermal expansion γ (3λ + 2μ)α i J n (.) The Bessel function of the first kind of order n u i Components of the displacement vector F i Component of the external forces per unit mass Q Intensity of the external heat sources per unit mass G i Component of the external applied couple per unit mass J i j Microrotation tensor μ ji Component of the moment of couple stresses ω j Vector of internal rotations q i Component of the heat flux vector e i jk Alternating tensor η Entropy

Introduction
Eringen [1], Nowacki [2], and Iesan [3] extended the equations of the micropolar theory of elasticity to include the thermal effects. Shanker and Dhaliwal solved several plane strain problems for an infinite body [4] and solved some dynamic thermoelastic problems in micropolar theory [5]. Chandrasekharaiah obtained the equations of heat-flux-dependent micropolar thermoelasticity, proved reciprocal principles and variational for his equations [6,7]. Sherief et al. [8] derived the generalized theory of micropolar thermoelasticity. This theory predicts a finite speed of propagation of the thermal and mechanical effects. Some contributions to this topic can be found in [9][10][11][12]. Many existing models of physical processes have been modified by fractional calculus. Caputo and Mainardi [13,14], Caputo [15], and Adolfsson et al. [16] used fractional derivatives to build good models of the behavior of viscoelastic materials.
In this manuscript, the theory of generalized micropolar thermoelasticity is modified. The equations of motion and constitutive relationships remain the same as the generalized theory of micropolar thermoelasticity. The modified energy equations are used to study two-dimensional axisymmetric temperature distribution problem in a half-space. Some comparisons are shown in Figures to estimate the impact of the fractional order parameter in all fields studied. All numerical results are in full agreement with previous work in the field of thermoelastic theories (see, for example, [8]).
The results presented in this research are critical in many scientific and engineering applications. The problems of cylindrical bodies are important because they have multiple applications in the production of mechanical parts.

Derivation of the governing equations
In this Section, we shall derive the governing equations of the fractional order theory of micropolar thermoelasticity. We shall use the Caputo fractional derivatives of order α ∈ [0, 1] of the continuous function f (t) given by [37]: where I β is the fractional integral of the function f (t) of order β defined by Miller and Ross [38] as: The equation of motion has the form [8] The balance equation of moments in the anisotropic case has the form [8] The second law of thermodynamics has the form The total energy of the system may be written as [8] where E 1 is the kinetic energy and E 2 is the internal energy of the system. We have where U is the internal energy per unit mass. The power of the external agents is given by From Eqs. (4)- (7), we obtain upon using Gauss's divergence theorem where Since the integrands in Eq. (8) are assumed to be continuous and the volume V is taken as arbitrary, Eq. (8) can be written in point wise form as: Interchanging the dummy indices i and j, we can write the preceding equation as Introducing Helmholtz's free energy function Eq. (11) takes the form Using the chain rule of differentiation, we can write Comparing Eqs. (13) and (14), we immediately obtain Expanding in Maclaurin's series up to quadric terms, we obtain In the preceding expression, the coefficient is all isothermal material constants for homogeneous materials, θ The following symmetries are evident: If we suppose that, in the reference configuration (θ 0), the material is free from stress and couple stress and has zero entropy, then from Eqs. (15)- (18) we can find the following restrictions on the material constants: Without loss of generality, we may also take 0 0. This is because the inclusion of 0 has no effect on the expressions (15)- (17).
ρη ρc E T 0 θ − a i j u i j − p i j ω i j (22) which are the constitutive equations of the theory of linear micropolar generalized thermoelasticity. To these constitutive equations, we add the fractional generalized Fourier's law of heat conduction of the form where τ 0 is a constant that acts as a relaxation time.
To obtain the equations of motion, we substitute from Eqs. (20) and (21) into Eqs. (1) and (2), and, upon Eqs. (9), (10), and (19), we get To obtain the equation of heat conduction in our case, we first write the linearized version of the equation of energy (3), which is By taking the time derivative of the both sides of the preceding equation, upon using Eq. (23) we obtain Substituting from Eqs. (22) into Eqs. (27), we obtain For the case of isotropic media, we have where ε m , α m , β m , and υ m are constants of the theory of micropolar thermoelasticity. Substituting from Eq. (29) into the constitutive Eqs. (20)-(23), these equations take the form Substituting from Eqs. (29.1)-(29.6) into Eqs. (24), (25), and (28), these equations become In obtaining Eq. (35), we have used the following relations: Substituting from Eqs. (9) and (10) for the quantities u ij and ω ij into Eqs. (34)-(36), we obtain In obtaining Eq. (38), we have used the fact that Equations (37) and (38) can be written in vector form as Substituting from Eqs. (9) and (10) into Eqs. (30) and (31), we obtain This completes the formulation of the fractional order theory of micropolar thermoelasticity. The governing equations consist of the equations of motion (40) and (41), and the equation of heat conduction (39). These equations are supplemented by the constitutive Eqs. (42) and (43).

An axisymmetric thermal shock problem for a half-space
In this Section, we shall obtain the solution for the problem of a half-space whose boundary is rigidly fixed and subjected to an axisymmetric thermal shock. There are no body forces, body couples, or heat sources affecting the medium.
We shall use a cylindrical coordinate system (r, ϕ, z) with the z-axis coinciding with the axis of symmetry. The surface of the half-space is taken as the plane z 0, with the z-axis pointing inward.
Due to symmetry, all the considered functions will be functions of the variables r and z only and will be independent of ϕ. Moreover, the components of the displacement and of the rotation vectors will have the form The cubical dilatation e is given by The operator ∇ 2 in the above equations is given by The governing equations have the form The constitutive equations have the form The boundary conditions of the problem state that the surface of this medium is rigidly fixed and is subjected to an axially symmetric thermal shock. The boundary conditions are thus The governing equations can be put in a more convenient form by using the non-dimensional variables where η 0 and v 1 are constants given by η 0 ρc E /k, v 1 √ (λ + 2μ)/ρ . Note that η 0 1/a where "a" is the thermal diffusivity coefficient.
The governing equations take the following forms where we have dropped the primes for convenience: The boundary conditions become ,

Solution in the Laplace transform domain
Taking the Laplace transform (denoted by an over bar) of both sides of Eqs. (58)-(60) and using the homogenous initial conditions, we get Making use of the relations Eqs. (71) and (72) take the forms Applying the divergence operator to both sides of Eq. (74), we get Eliminating e between Eqs. (73) and (76), we obtain In a similar manner we can show that e satisfies Eq. (77), i.e., Equation (77) can be written in a factorized form as: where k 2 1 and k 2 2 are the roots with positive real parts of the characteristic equation The solution of Eq. (78) can be written as where We shall use the Hankel transform defined by the relation [39,40] where J 0 is the Bessel function of the first kind of order zero. This transform has the inversion formula [39,40] f (r , z, s) The operational property of the Hankel transform is [39,40] Applying the Hankel transform to both sides of Eq. (81), we obtain where the operator denotes partial differentiation with respect to z. Thus, the solution of Eq. (82) has the form where A 1 and A 2 are parameters that depend on q, s, and q i is given by Similarly, the solution of Eq. (77) compatible with Eq. (75) has the form: Applying the inversion formula of the Hankel transform to Eqs. (83) and (84), we get The components of Eqs. (73) and (74) are: Differentiating Eq. (44) with respect to r and substituting the result into Eq. (87), we get Similarly, we can write Eq. (88) in the form Letting Eqs. (87), (88), and (89) take the forms: Differentiating Eq. (93) with respect to z and Eq. (94) with respect to r the difference between the resulting equations takes the form: We now take: Substituting from Eqs. (97) into Eqs. (95) and (96), and using the fact that we get upon integrating with respect to r Eliminating H between Eqs. (98) and (99), we obtain In a similar manner we can show that H satisfies the equation Equation (100) can be factorized as where n 2 1 and n 2 2 are the roots of the equation As before, we have the solution * where D 1 and D 2 are parameters that depend on q, s, and m i is given by The inverse formula of the Hankel transform of Eq. (104) after using Eqs. (97) is where we have used the relation [40] d dr J 0 (qr) −q J 1 (qr).
Similarly, the solution of Eq. (101) gives Hankel-Laplace transform of H , in the form The inverse formula of the Hankel transform of Eq. (106) after using Eq. (97) is In obtaining the preceding expressions, we have used the following relation [41,42]: We shall apply the boundary conditions of the problem to find the parameters A 1 , A 2 , D 1 , and D 2 . The transformed boundary conditions (57) have the form: From boundary conditions (118), (119) and Eqs. (85), (105), (108), and (109), respectively, we get Solving the linear system of Eqs. (120)-(123), numerically, we obtain the parameters A 1 , A 2 , D 1 , and D 2 . This completes the solution in the Laplace transform domain.

Numerical results
The polystyrene (one of the polymers) was chosen for purposes of numerical evaluations. The constants of the problem were taken as in Table 1.
To compute the values of the functions, a numerical procedure was used to invert the double transforms in the above expressions. First, a numerical method based on Fourier expansion was used to invert the Laplace transform [43]. Next, the integrals appearing in the formula for the inverse of the Hankel transform were evaluated numerically by a subroutine which is a modification of "qint" from the book [44]. The FORTRAN programming language was used throughout on a personal computer. The accuracy maintained was five digits.
The non-dimensional temperature of the surface is taken as where θ 0 is a constant and H is the Heaviside unit step function. This means that starting at time t 0, a circle of unit radius is suddenly raised to a temperature equal to 00 and kept at this temperature, while the rest of the surface is kept at zero temperature. Taking the Hankel transforms, we obtain The constant θ 0 is taken equal to unity during computations.  To study different aspects of the solution, we did three groups of graphs. Figures 1, 2, 3 show the solution on the z-axis (r 0) as function of z for t 0.1 and different values of α {0, 0.9, 1}. Figures 4, 5, 6 show the solution on the z-axis also as function of z but for a fixed α 0.8 and different values of t {0.05, 0. 1}. Figures 7,8,9,10,11,12,13 show the solution on the plane z 0.15 as function of r for a fixed α 0.5 and different values of t {0.5, 1.0}. Figures 1, 2, 3 show the difference between coupled, generalized, and fractional thermoelastic micropolar theories. For α 0 and τ 0 0, the solution is that of the coupled theory where the speed of the wave is infinite but for α 1, the solution is that of the generalized theory where the speed of the wave is finite. For α close to 1, the solution seems to behave like that of the generalized theory.
For the generalized case, we have four waves. The fastest wave is mainly thermal. Its wave front in Figs. 1, 2, 3 has reached the point z 0.71, approximately, meaning a wave speed of 7.1 [8]. Outside the region 0 < z < 0.71, all functions are zero. This is not the case for the coupled and uncoupled theories where an infinite For the fractional case, the fast wave front has reached the point z 0.861 meaning a wave speed of 8.61. Outside the region 0 < z < 0.861, all functions seem to be zero within the accuracy maintained during computation. This is an indication that the wave has finite speed of propagation, not a proof. The authors are working on a theoretical proof of this fact along the lines in [31]. The second wave front occurs also at z 0.1 approximately.

Conclusions
We have derived a mathematical model to describe the interaction of mechanical and thermal effects on a micropolar thermoelastic medium using the methodology of fractional calculus. The main advantage of the fractional calculus model over that of the generalized model derived by Sherief et al. [31] is that it predicts a retarded response to these effects as is the case in nature as opposed to the instantaneous response of the generalized theory.
Povstenko [45] reports that in the case of the time fractional diffusion equation with (in our notations) 0 < α < 0.36 the solution has features more closely resembling those of the heat conduction equation. For 0.69 < α < 1, the solution resembles the solution to the wave equation, and the intermediate behavior corresponds to the values 0.36 < α < 0.69 (see also Fujita [46]).
Of course, the situation in the case of fractional thermoelasticity is much more complicated. The numbers in the last paragraph will be functions of the material parameters of the medium. Many numerical experiments were performed. They indicate, roughly speaking, that for α close to zero, the solution behaves similar to that of the coupled theory. For α close to unity, the solution behaves similar to that of the generalized theory.
We like to point out here that the behavior of the solution depends on α and also on the type of fractional derivative used. All the above statements are valid for the original Caputo fractional derivative. In [31], Sherief and Raslan proved that for a 1D problem when using the fractional theory with the new Caputo definition of fractional derivative then we have two waves. The slower wave moves with finite speed while the other wave moves with infinite speed for all values of α < 1. This means that the solution resembles that of the coupled theory for all values of α < 1.
This manuscript introduces, for the first time, a new direct approach to solve 2D problems of micropolar thermoelasticity (valid for both the generalized and the fractional theories). This avoids the use of potential functions with their known disadvantages.
We see from Figs. 1, 2, 3 that a change in α has a small effect on the temperature and a larger effect on the displacement and stress components. Figures 1, 2, 3, 4, 5, 6 show that the temperature and stress have discontinuities at the wave fronts but the displacement, as expected, is always continuous. We know theoretically that there are four waves. The effect of the other two waves is somewhat small and cannot be studied numerically. The authors are working on a method to study these waves theoretically as was done in [31].

Data availability statement
The manuscript has no associated data.

Conflict of interest There is no conflict.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.